初心者向きにR言語を用いたIIRフィルターの作成方法をまとめました。記事が冗長ですが、お許しください。
IIRフィルターは、ノイズ除去や特定の周波数帯域の信号抽出でよく用いられます。ローパスフィルターやバンドパスフィルター等として使うことができます。今回はR言語を用いて実装しました。
パッケージの準備
今回は signal パッケージ[1]を使用し、IIRフィルターを作成します。まず、パッケージをインストールします。
#signalパッケージのダウンロード
install.packages('signal')
#ライブラリーの読み込み
library('signal')
入力データは10Hzと50Hzの正弦波の合成波を使用します。
#入力データの作成（周波数が10Hzと50Hzのサイン波のデータの合成、サンプリング周波数：1kHz、１０秒間のデータ）
time <- seq(from=0, by=1/1000, length=10000) #時間軸
data <- sin(2*pi*10*time) + sin(2*pi*50*time)
#入力データの図示
plot(time, data, t="l", xlim=c(1, 1.5))
上図より、2つの周波数の合成波になっています。また、周波数解析を行うと、
data_fft <- fft(data)
data_psd <- abs(data_fft)^2
plot_fft <- seq(from=1/10000, by=1000/10000, to=500)
#入力データの周波数特性の図示
plot(plot_fft, data_psd[1:(length(data_psd)/2)], t="l", log="x", xlab="frequency (Hz)", ylab="PSD (1/Hz)")
ピークが2つ（10Hzと50Hz）になっていることがわかります。
IIRフィルターの設計
今回はIIRフィルターの一種であるButterworth filterを実装するために、signalパッケージのbutter()を使用します。帯域通過周波数が5〜15Hzのバンドパスフィルターを実装します。
#バンドパスフィルタの設計
#バターワースフィルターの各パラメータ
fil_N <- 2 #フィルターの次数
fs <- 1000 #サンプリング周波数
fn <- fs/2 #ナイキスト周波数
fc <- c(5,15) #帯域通過周波数
#fc <- 5 #ローパスフィルターやハイパスフィルターでは、周波数の指定は一つの値のみ
fc_norm <- fc/fn #周波数の正規化
butterworth_filter <- butter(fil_N, fc_norm, type="pass", plane="z")
butter()では、帯域周波数をナイキスト周波数で割った正規化周波数を使います。typeは、pass, low, high, stopの指定により、様々な種類のフィルターを実装できます。
フィルター特性の確認
freqz(butterworth_filter, Fs=fs) #フィルター特性の図示
5〜15Hzのバンドパスフィルターを設計できました。
設計したフィルターを入力データに通す
#設計したフィルターを入力データに適応
data_fil <- filtfilt(butterworth_filter, data)
フィルターによる位相遅れを防ぐためにfiltfilt()を使用し、順方向と逆方向の両方の処理を行うゼロ位相デジタルフィルター処理を実装します。
フィルターを通したデータの確認
#フィルターを通したデータの図示
plot(time, data_fil, t="l", xlim=c(1, 1.5))
#周波数特性の図示
data_fft <- fft(data_fil)
data_psd <- abs(data_fft)^2
plot(plot_fft, data_psd[1:(length(data_psd)/2)], t="l", log="x", xlab="frequency (Hz)", ylab="PSD (1/Hz)")
フィルターにより50Hz成分がなくなりました。
参考文献
[1] Uwe Ligges et al. "Signal Processing", ver. 0.7-6, https://cran.r-project.org/web/packages/signal/signal.pdf, 参照Sept 21, 2018.