scilab_FIRフィルタ

clear
// ファイル読み込み
[s,Fs,bits]=wavread('sample.wav');
// wavのサイズ
N = size(s,'c');

// ハミング窓関数生成
win = window('hm',N);

// フーリエ変換
FFT = fft(s.*win);
// 振幅特性を測るため絶対値をとる
FFT_abs = abs(FFT);

// x軸の設定
// Fs/N：分解能、振幅スペクトルはN/2で線対称となるので0~N/2の範囲を使う
x = Fs*(0:(N/2))/N;

// グラフ出力
scf(0);
clf;
subplot(2,1,1);
plot(s);
scf(1);
clf;
subplot(2,1,1);
//縦軸そのまま
plot(x,FFT_abs(1:size(x,'c')));
xtitle('FFT分析','[Hz]','');

//フィルタ生成
[h,hm,fr]=wfir('sb',701,[540/Fs 650/Fs],'hm',[1 0]); //ハミング窓

//フィルタの適用(畳み込み)
snd=convol(s,h);

scf(0);
subplot(2,1,2);
plot(snd);

scf(2);
clf;
// システム関数
plot(fr,hm); // 正規化周波数で横軸をとる
scf(3);
clf;
// インパルス応答
plot(h);

// wavのサイズ
Naf = size(snd,'c');

// ハミング窓関数生成
winaf = window('hm',Naf);

// フーリエ変換
FFTaf = fft(snd.*winaf);
// 振幅特性を測るため絶対値をとる
FFTaf_abs = abs(FFTaf);

// x軸の設定
// Fs/N：分解能、振幅スペクトルはN/2で線対称となるので0~N/2の範囲を使う
xaf = Fs*(0:(Naf/2))/Naf;

scf(1);
subplot(2,1,2);
plot(xaf,FFTaf_abs(1:size(xaf,'c')));
xtitle('','[Hz]','');

//サウンド保存
savewave('result.wav',snd);

FIRフィルタをscilabで設計した。音声ファイルを読み込んで解析した後真ん中の周波数をカットするもの。