A13 FIRフィルタ その2
室長：測定編のまえに、FIRフィルタの基本的な事をおさえておくぞ。いくらこの研究室がなるべく理論に触れないといっても最低限必要な知識があるからな。
助手：難しい数学はいやですよ。
室長：大丈夫だ、問題ない。簡単な事しかいわん。まず入力と出力が全く同じになるフィルタはどんなフィルタかな？
助手：これは簡単です。係数が掛け算されて足されるのだから係数に１が１つだけあるフィルタです。たとえば[ 1 0 0 0 ]とかです。
室長：うむ。またフィルタの前にゼロをつけても出力は同じというのはわかるな。[ 0 1 0 0 0 ] 出力される時間がゼロの数ぶん遅れる。また後ろにゼロをつけても出力結果は変わらない。あとフィルタの出力を半分にするにはフィルタの係数に0.5をかければ良い。
助手：わかります。
室長：あと直線位相のフィルタの形は左右対称なのじゃ。理論的な事は他のサイトをあたってくれ。
この左右対称な直線位相フィルタの前にゼロを挿入して
こうしても時間遅れがあるだけでやっぱり直線位相フィルタじゃ。
助手：なるほど、係数のなかの位置が変わっても波形自体が左右非対象なら直線位相なのですね。
室長：あとフィルタの次数が偶数か奇数かで特性が変わってくる。フィルタの次数が偶数（tap数は奇数）でないと左右対称の直線位相ハイパスフィルタが作れないのじゃ。なぜならtap数が偶数で左右対称にしようとすると必ずピークの幅が2サンプルになってしまうからじゃ。
fs=9600;
f=[0 300 310 fs/2];
m=[-300 -300 0 0];
f=f/(fs/2);
m=db2a(m);
tap=512;
a=fir2(tap-1,f,m,8192,kaiser(tap,8));
tap=511;
b=fir2(tap-1,f,m,8192,kaiser(tap,8));
b=[b;0]; #bの後ろのゼロを1つ追加
clf;
hold on;
plot(a);
plot(b,'1');
助手：tap数が1違うだけでずいぶん波形が変わりますね。
室長：この波形の変化を嫌うならtap数を奇数で設計して、後ろにゼロを１つ付ければよい。brutefirのtap数は2の乗数と決まっているからな。1/2サンプル分出力時間がずれるのでやるなら全部のフィルタでやった方が良いじゃろう。音を聴いて違いがわかるかは試していないがの。
あと信号にフィルタをかけたときの波形の見方をやってみよう。octaveのfftconv関数を使う。
助手：はいはい。いつも通り’help fftconv’をみてと。
a=[zeros(100,1) ; ones(441,1) ; -ones(441,1) ; zeros(100,1)];
b=fir1(510, 500/22050, "low");
c=fir1(510, 500/22050, "high");
ab=fftconv(a,b);
ac=fftconv(a,c);
clf;
hold on;
plot(a,'1');
plot(ab,'2');
plot(ac,'3');
1kHzの矩形波の波形1周期分を500Hzのローパスフィルタ、ハイパスフィルタにかけてみました。
ところでfftconvではaとbの順序がどちらでも良いようですが？
室長：実はな、フィルタの係数とデータの波形は同じ扱いなのじゃ。例えばこれから測定するインパルスレスポンスはそのままフィルタの係数として使える。またいままで使ってきた自作のfilterplot関数はフィルタ波形をfftして周波数成分を表示していたわけじゃ。
じゃあ500Hzのハイパスフィルタと2kHzのローパスフィルタを畳み込むとどうなるかな。
助手：順序はどうでもよいですが、ハイパスフィルタを通した後にローパスフィルタを通すイメージですね。ということは500Hzから２kHzのバンドパスフィルタになるはずです。
a=fir1(510, 2000/22050, "low");
b=fir1(510, 500/22050, "high");
c=fftconv(a,b);
clf;
subplot(2,1,1);
plot(c);
subplot(2,1,2);
filterplot(c,44100);
室長：あらかじめわかっているのなら最初からバンドパスフィルタを設計したほうが面倒が無いがのう。 さて次回だが測定の前に一回寄り道じゃ。
2012年12月17日