FFTコード（説明）

えー先日勢いでアップしましたFFTコードですが
やっぱり説明があったほうが親切なカンジが出て好感触なんじゃないかと
今思いましたので、説明を入れてみようかなと。

■FFT実行ライブラリ go_fft
入力された信号に対してfft結果を返すライブラリです。

---

function [amp, frequency]=go_fft(input_data, fs, win_sel, out_form)

*引数
・input_data : 入力信号列の変数の指定。
・fs : 入力信号のサンプリングレートの指定。
        ちなみに'1'を指定すれば結果的に正規化周波数でのFFT実行になります。
・win_sel : 窓関数の指定。対応している窓関数は以下の３つのみです。
             're' : 矩形窓。窓関数を使用しない場合に指定ください。
             'hn' : hanning窓を指定
             'hm' : hamming窓を指定
・out_form : 出力フォーマットを指定。
              'db' : FFT結果の値"1.0"を基準値としたdB単位での出力。
                      つまり、ある周波数でのFFT結果が1.0であった場合'0[dB]', 
                   0.5であったら'-6[dB]'の値が出力されます。
              'lnr' : FFT結果をリニア値で出力します。
                       上記の指定'db'のｄＢ変換する前の値となります。
              'abs' : FFT結果の生データである複素数値を絶対値に変換した値を返します。
                       上記の'db', 'lnr'の値は、窓関数の補正を施した値ですが、
                       'absの場合、その補正が入っていない値となります。
              'cmp' : fft結果の生データである複素数を返します。
                       具体的には、scilabの関数fft()の出力値となります。 
*戻り値
・amp : FFT結果である各周波数ごとの振幅値の配列です。
         この値の形式は引数'out_form'に依存します。
・freq : FFT結果である各binの周波数を示す配列です。
          この配列'freq'の並びは'amp'の並びに対応しています。 

---

■go_fft実行コード

//-----------------------------------
// sample code to execute go_fft
// ver: 1.0 
//-----------------------------------

//set lib path
//Please change the following path to fit your enviroment. 
getf('C:\work\scilab\go_fft.sci');  

ここでは、ライブラリ'go_fft'のファイルをリンクしています。
みなさまの環境に合わせて、'go_fft,sci'が格納されている
ディレクトリにご変更ください。

//generate test wave
fs = 10000;
f_1 = 1250;
f_2 = 2500;
n = [1:1:8192];
sine_1 = 0.5 * sin(2*%pi*(f_1/fs)*n);
sine_2 = 0.25 * sin(2*%pi*(f_2/fs)*n);
input_wave = sine_1 + sine_2; 

ここでは、FFTをかける信号を生成しております。
1250[Hz]、振幅0.5のsine波と
2500[Hz]、振幅0.25のsine波とを加算したものを基準信号としています。

//go fft
[amp f] = go_fft(input_wave, fs, 're', 'lnr'); //plot in linear unit
//[amp f] = go_fft(input_wave, fs, 're', 'dB'); //plot in dB unit 
//[amp f] = go_fft(input_wave, fs, 'hn', 'lnr'); //use hanning window
//[amp f] = go_fft(input_wave, fs, 'hm', 'lnr'); //use hamming window 

ライブラリ'go_fft'のコール、窓関数は矩形（窓関数なし）、でdBフォーマットを指定しています。

//display fft result
plot(f, amp);
xgrid; //grid on 

結果のグラフ表示です。

---

■go_fft本体

//---------------------------------------
// go_fft.sci
// calculate fft on scilab
// ver : 1.0
//--------------------------------------- 

function [amp, frequency]=go_fft(input_data, fs, win_sel, out_form) 

ライブラリの引数戻り値指定。=の隣にスペースを入れるとコンパイルできない？ようで、
数時間カッツリはまったストレスからウィスキーを痛飲したことは、今ではいい思い出です。

printf("===== go_fft =====\n");
printf(" window   = %s\n", win_sel);
printf(" out form = %s\n", out_form); 

//get input length
input_length = length(input_data); 

printf(" input_length = %d\n", input_length);

入力信号配列の長さを取得。

//output length
out_length = input_length/2; 

出力信号の長さを決定。折り返し周波数のfs/2以上の周波数は出力しません。

//windowing
if win_sel == 'hn'
  window_adjust_coef = 0.5;
elseif win_sel == 'hm'
  window_adjust_coef = 0.54;
else
  window_adjust_coef = 1;
end 

窓関数を使用した場合、その影響でFFT結果が小さくなるため
補正値をFFT結果にかける必要があります。その補正値の指定ッス。

win_l = window(win_sel, input_length); 

printf(" window_adjust_coef = %f\n",window_adjust_coef);
  
if win_sel ~= 're'  
  input_windowing = input_data .* win_l;
else
  input_windowing = input_data;
end 

窓関数を入力にかけてます。

//fft process
res_cmp = fft(input_windowing);

FFTの実行

//res_cmp = fft(input_data);
res_abs = abs(res_cmp)*(2/input_length);

FFT結果である複素数を絶対値に変換

res_lnr = res_abs/window_adjust_coef; 

窓関数の補正実行

//convert dB unit
res_db = 20 * log10(res_lnr); 

FFT結果のdB単位への変換

if out_form == 'dB'  
  result = res_db;
elseif out_form == 'lnr'
  result = res_lnr;
elseif out_form == 'abs'
  result = res_abs;
elseif out_form == 'cmp'
  result = res_cmp;
else
  result = res_db;
end

amp = result(1:out_length);

引数で指定された形式の値を戻り値'amp'に代入

//generate frequency 
freq_points = [0:1:out_length-1];
frequency = (fs/input_length) * freq_points;  

戻り値'freq'を'fs'と入力信号の長さから生成

printf("=================="); 

endfunction 

---

以上ッス。