前回、FIRフィルタを作ろうとしたら、いくつも設計用の関数が並立？していることに気づきました。その中でfsfirlin関数というものの設計例通りに「手習ひ」しバンドパスフィルタらしきものを生成しました。しかし、実際に信号をフィルタしてません。信号処理素人の老人は実際の波形に適用してみないと納得いかんね。
※「手習ひデジタル信号処理」投稿順 Indexはこちら
※Windows11上で、Scilab6.1.1およびScilab上のツールボックス Scilab Communication Toolbox 0.3.1（以下comm_tbx)を使用させていただいとります。
以前に定義済の関数再掲載
計算してみるのに過去回で作った自前関数を２つばかり使ってます。いつ頃作ったかも忘却力の老人には定かでないので、再掲載しておくことにいたしました。最初はサンプル波形用にSIN波を定義するためのもの。
// Generate Sine wave
// Fs: sampling frequency [Hz]
// nSample: number of samples
// fq: frequency [Hz]
// phdeg: phase [deg]
// mag: manitude
function [t,y]=genSin(Fs, nSample, fq, phdeg, mag)
t = ((1:nSample) - 1) / Fs;
y = mag*sin(2*%pi*(fq*t+(phdeg/360)))
endfunction
つづいてFFTプロット（縦軸対数）を行うもの
// FFT plot (dB)
// Fs: sampling frequency [Hz]
// wav: time domain seq(nSample)
// opt: Window function 0:none, 1:hamming
function plotFFT2(Fs, wav, opt)
N = size(wav, '*');
regN = 2/N;
select opt
case 1 then
winf = window('hm',N);
regF = 2;
else
winf = ones(1,N);
regF = 1;
end
y = fft(wav .* winf);
xf = Fs*(0:(N/2))/N; //associated frequency vector
nf = size(xf, '*');
clf();
plot2d(xf, 20*log10(abs(y(1:nf)) * regN * regF));
xtitle('FFT', '[Hz]', '[dB]');
endfunction
実験してみるフィルタの定義
fsfirlin関数は「こういう特性にしてね」と標本（目標）ベクトルで与えてお願いすると、FIRフィルタを設計してくれる関数のようです。そこで今回実験してみるのは以下のような入力定義から生成できるローパスフィルタ（LPF）です。
-
- fsfirlinの入力となる周波数応答標本のサンプル数は６４点（前回と同じ）
- 標本「８」までは通過域
- 標本「９」で振幅半減（よってカットオフ周波数は標本「８」と「９」の中間、正規化周波数で0.13付近にあるはず）
- 標本「１０」以降は０
- ２型
なお、与えるサンプル点数を長くすると、生成されるフィルタの次数も長くなるみたいデス。知らんけど。上記の場合のフィルタ定義は以下に。
hd=[ones(1,8) zeros(1,56)];
hd(9)=.5;
hst=fsfirlin(hd,2);
pas=1/prod(size(hst))*.5;
fg=0:pas:.5;
plot2d(fg(1:257)',hst);
xtitle('LPF(fsfirlin)', 'Normarized Frequency', 'Magnitude');
上記で生成されたフィルタの特性のプロット。
実験材料の波形定義
通過域と阻止域にサイン波を一つづつおいたものを入力波形として、これに上記のフィルタを適用してみることにいたしました。入力波形は以下のようです。
-
- サンプリング周波数 10kHz
- サンプル点数4096
- 第１周波数（通過するハズ）440Hz
- 第２周波数（阻止されるハズ）3333Hz
Fs=10e3; nSample=4096; H440=440; H3333=3333; [t, y440]=genSin(Fs, nSample, H440, 0, 1); [t, y3333]=genSin(Fs, nSample, H3333, 0, 1); wavT = y440 + y3333; clf; plot2d3(t, wavT); clf; plotFFT2(Fs, wavT, 1);
生成された時間波形の一部を拡大したものが以下に。
上記波形のFFTプロットが以下に。
フィルタしてみる
実験してみたコードが以下に。こんなんで良いのか？
y=filter(hst, 1, wavT);
plot2d3(t, y);
clf
plotFFT2(Fs, y, 1);
xtitle('Output(LPF:fsfirlin)', 'Frequency[Hz]', 'Mag[dB]');
まずは時間波形。
そしてFFTの結果。
いいんかな～、こんなことで。闇雲。