RでFIRフィルタ設計(2)

Rでfir1を使って信号をフィルタリングしてみる
R関連の情報は少ないので、MATLAB関連のページを見て、実際にRに適用してみる。

http://www.phon.ucl.ac.uk/courses/spsci/matlab/lect8.html

上記ページを参考に、Rで10000Hzでホワイトノイズを作成して、2000-3000Hzをフィルタリングするバンドパスフィルタを作成する。

フィルタ対象のホワイトノイズを作成する

> Fs = 10000 # サンプリング周波数
> x = rnorm(10000)
> plot (x)

バンドパスフィルタを作成する

> N = 128 #フィルタの次数
> fil = fir1(N-1,c(2000,3000) / (Fs / 2), type="pass")

第２引数のc(0,1)の1はナイキスト周波数（Fs/2）となるので、Fs/2でスケーリングする。

このフィルタの周波数特性を確認する

> freqz(fil, Fs=Fs)

2000-3000Hzのところにフィルタができていることを確認できる。

ホワイトノイズを時間領域で畳み込み、フィルタリングする

> x2 = convolve(x,fil,type="o")

ここでは、時間領域の畳み込みのために、convolve()を使う。conv()でもいいみたいだが、conv()の戻り値はtsとして返ってくるので扱いに注意が必要。
convolve()のソースを見ると入力信号の２つをフーリエ変換して、積をとっている。
時間領域の畳み込みは周波数領域での積となり、周波数領域ではその逆となる（らしい）。
http://www.wakayama-u.ac.jp/~kawahara/signalproc/CH2cont2discrete/signalshint.html#anchor1831284
ちなみに、typeは”open”にしておく。"circular"や”filter”では所望の信号を得ることができない。

あと、畳み込みにより、信号長が変わっている。これは畳み込みで、n = length(x) + length(fil) - 1の長さで処理が行われるためである。

> length(x)
[1] 10000
> length(x2)
[1] 10127

なので、信号を128から10127の区間にする。

> x3 = x2[128:length(x2)]

フィルタリング結果を確認する

> plot(x3)

これだとあまりよくわからないが、振幅に変化があることはわかる。

フィルタリング結果のパワースペクトラムを見てみる

> plot(abs(fft(x3))^2,type="l")

2000-300Hzでフィルタリングされていることがわかる。
7000-800Hzのところもスペクトルがあるのは、サンプリング定理により、ナイキスト周波数(5000(=Fs/2))でエイリアシングがおきているため。

最後に

なんとか、Rでバンドパスフィルタの作成とその適用をすることができた（ように思う）。あまりにも、これ系のページが少ない＋難しくてわかりにくいので、わかりやすいようにやりかたを書いてみた。間違いがあるかもしれないので、そのときは、コメント欄にでも指摘よろしく。