ローパス FIR フィルター – ウィンドウ法

この例では、FIR フィルターを設計し実装する方法を示します。2 つのコマンド ライン関数fir1とdesignfiltおよび対話型のフィルターデザイナーアプリを使用します。

例で使用する信号を作成します。信号は、加法性 $$N(0,1/4)$$ホワイト ガウス ノイズを伴う 100 Hz 正弦波です。再現性のある結果を得るために乱数発生器を既定の状態に設定します。

rngdefaultFs = 1000; t = linspace(0,1,Fs); x = cos(2*pi*100*t)+0.5*randn(size(t));

フィルター設計は、次数が 20、カットオフ周波数が 150 Hz の FIR ローパス フィルターです。長さがフィルター次数よりも 1 サンプル大きく、 $$\beta =3$$のカイザー ウィンドウを使用します。カイザー ウィンドウについての詳細は、関数kaiserを参照してください。

fir1を使用してフィルターを設計します。fir1は、区間 [0,1] で正規化された周波数を必要とします。ここで、1 は $$\pi$$ラジアン/サンプルに相当します。fir1を使用するには、すべての周波数仕様を正規化周波数に変換しなければなりません。

フィルターを設計してそのフィルターの振幅応答を表示します。

fc = 150; Wn = (2/Fs)*fc; b = fir1(20,Wn,'low',kaiser(21,3)); fvtool(b,1,'Fs',Fs)

信号にフィルターを適用し、100 Hz の正弦波の最初の 10 周期の結果をプロットします。

y = filter(b,1,x); plot(t,x,t,y) xlim([0 0.1]) xlabel('Time (s)') ylabel('Amplitude') legend('Original Signal','Filtered Data')

designfiltを使用して同じフィルターを設計します。フィルター応答を'lowpassfir'に設定し、仕様をName,Valueのペアとして入力します。designfiltでは、フィルターの設計を Hz で指定できます。

Fs = 1000; Hd = designfilt('lowpassfir','FilterOrder',20,'CutoffFrequency',150,...'DesignMethod','window','Window',{@kaiser,3},'SampleRate',Fs);

データをフィルター処理して結果をプロットします。

y1 = filter(Hd,x); plot(t,x,t,y1) xlim([0 0.1]) xlabel('Time (s)') ylabel('Amplitude') legend('Original Signal','Filtered Data')

フィルターデザイナーを使用したローパス FIR フィルター

この例では、ウィンドウ法を対話型のフィルターデザイナーアプリで使用して、ローパス FIR フィルターの設計と実装を行う方法を示します。

y2 = filter(Hd,x); plot(t,x,t,y2) xlim([0 0.1]) xlabel('Time (s)') ylabel('Amplitude') legend('Original Signal','Filtered Data')

また、対話型のツールfilterBuilderを使用してフィルターを設計することもできます。

バンドパス フィルター – 最小次数の FIR および IIR システム

この例では、バンドパス フィルターを設計して、最小次数の FIR 等リップル フィルターと IIR バタワース フィルターを使ってデータをフィルターする方法を示します。実際の信号の多くは、振動コンポーネント、低周波のトレンドおよび加法性ノイズの重ね合わせとしてモデル化できます。たとえば、経済データには多くの場合、ゆっくり変化する上向きまたは下向きのトレンド上に重ね合わされた周期を表す変動が含まれています。さらに加法性ノイズ コンポーネントにより、変動過程における測定誤差と固有のランダムな変動の組み合わせが加えられます。

以下の例では、何らかの過程を 1 年間毎日サンプリングすると仮定します。過程は、概ね週単位、月単位での変動すると仮定します。さらに、データは低周波の上向きトレンドをもち、加法性 $$N(0,1/4)$$ホワイト ガウス ノイズがあります。

信号を、1/7 および 1/30 周期/日の周波数をもつ 2 つの正弦波の重ね合わせとして作成します。低周波の増加トレンド項と $$N(0,1/4)$$ホワイト ガウス ノイズを追加します。再現可能な結果が必要な場合は、乱数発生器をリセットします。データは 1 サンプル/日でサンプリングされます。結果の信号とパワー スペクトル密度 (PSD) 推定をプロットします。

rngdefaultFs = 1; n = 1:365; x = cos(2*pi*(1/7)*n)+cos(2*pi*(1/30)*n-pi/4); trend = 3*sin(2*pi*(1/1480)*n); y = x+trend+0.5*randn(size(n)); [pxx,f] = periodogram(y,[],[],Fs); subplot(2,1,1) plot(n,y) xlim([1 365]) xlabel('Days') grid subplot(2,1,2) plot(f,10*log10(pxx)) xlabel(“周期/天”) ylabel('dB') grid

低周波のトレンドが、低周波パワーの増加としてパワー スペクトル密度推定に現れます。低周波パワーは、1/30 周期/日において変動の約 10 dB 上に現れます。この情報をフィルターの阻止帯域の仕様に使用します。

最小次数の FIR 等リップル フィルターと IIR バタワース フィルターを次の仕様で設計します。通過帯域が [1/40,1/4] 周期/日、阻止帯域が [0,1/60] および [1/4,1/2] 周期/日とします。両方の阻止帯域の減衰量を 10 dB に、通過帯域リップルの許容誤差を 1 dB に設定します。

Hd1 = designfilt('bandpassfir',...'StopbandFrequency1',1/60,'PassbandFrequency1',1/40,...'PassbandFrequency2',1/4 ,'StopbandFrequency2',1/2 ,...'StopbandAttenuation1',10,'PassbandRipple',1,...'StopbandAttenuation2',10,'DesignMethod','equiripple','SampleRate',Fs); Hd2 = designfilt('bandpassiir',...'StopbandFrequency1',1/60,'PassbandFrequency1',1/40,...'PassbandFrequency2',1/4 ,'StopbandFrequency2',1/2 ,...'StopbandAttenuation1',10,'PassbandRipple',1,...'StopbandAttenuation2',10,'DesignMethod','butter','SampleRate',Fs);

FIR および IIR フィルターの次数とアンラップされた位相応答を比較します。

fprintf('The order of the FIR filter is %d\n',filtord(Hd1))

The order of the FIR filter is 78

fprintf('The order of the IIR filter is %d\n',filtord(Hd2))

The order of the IIR filter is 8

[phifir,w] = phasez(Hd1,[],1); [phiiir,w] = phasez(Hd2,[],1); figure plot(w,unwrap(phifir)) holdonplot(w,unwrap(phiiir)) holdoffxlabel('Cycles/Day') ylabel('Radians') legend('FIR Equiripple Filter','IIR Butterworth Filter') grid

IIR フィルターは FIR フィルターよりもずっと低い次数をもっています。しかし、FIR フィルターには通過帯域にわたっての線形位相応答があり、これは IIR フィルターにはありません。FIR フィルターはフィルター通過帯域内のすべての周波数を同等に遅延しますが、IIR フィルターはそうではありません。

また、周波数の単位あたりの位相の変化率は、IIR フィルターよりも FIR フィルターでの方が大きくなっています。

比較用に、ローパス FIR 等リップル フィルターを設計します。ローパス フィルターの仕様は、通過帯域を [0,1/4] 周期/日、阻止帯域の減衰量を 10 dB 相当、通過帯域リップル許容誤差を 1 dB に設定するものとします。

Hdlow = designfilt('lowpassfir',...'PassbandFrequency',1/4,'StopbandFrequency',1/2,...'PassbandRipple',1,'StopbandAttenuation',10,...'DesignMethod','equiripple','SampleRate',1);

バンドパス フィルターおよびローパス フィルターを使ってデータをフィルター処理します。

yfir = filter(Hd1,y); yiir = filter(Hd2,y); ylow = filter(Hdlow,y);

バンドパス IIR フィルター出力の PSD 推定をプロットします。次のコードでyiirをyfirと置き換えて、FIR バンドパス フィルター出力の PSD 推定を表示できます。

[pxx,f] = periodogram(yiir,[],[],Fs); plot(f,10*log10(pxx)) xlabel(“周期/天”) ylabel('dB') grid

PSD 推定は、バンドパス フィルターにより低周波トレンドと高周波ノイズが減衰していることを示しています。

FIR フィルターと IIR フィルターの最初の 120 日間の出力をプロットします。

plot(n,yfir,n,yiir) axis([1 120 -2.8 2.8]) xlabel('Days') legend('FIR bandpass filter output','IIR bandpass filter output',...'Location','SouthEast')

FIR フィルターでの位相遅延の増大は、フィルター出力で明確に示されています。

比較のために、7 日および 30 日周期の重ね合わせの上にローパス FIR フィルター出力を重ねてプロットします。

plot(n,x,n,ylow) xlim([1 365]) xlabel('Days') legend('7-day and 30-day cycles','FIR lowpass filter output',...'Location','NorthWest')

上記のプロットでは、ローパス フィルター出力で低周波トレンドが明確なことがわかります。ローパス フィルターでは 7 日および 30 日の周期が保持されますが、バンドパス フィルターは低周波トレンドも削除するため、この例ではバンドパス フィルターの方が優れた挙動を示しています。

ゼロ位相フィルター処理

この例では、ゼロ位相フィルター処理を実行する方法を示します。

fir1とdesignfiltを使って、信号の生成とローパス フィルターの設計を再度行います。これらの変数が既にワークプレースにある場合は、次のコードを実行する必要はありません。

rngdefaultFs = 1000; t = linspace(0,1,Fs); x = cos(2*pi*100*t)+0.5*randn(size(t));% Using fir1fc = 150; Wn = (2/Fs)*fc; b = fir1(20,Wn,'low',kaiser(21,3));% Using designfiltHd = designfilt('lowpassfir','FilterOrder',20,'CutoffFrequency',150,...'DesignMethod','window','Window',{@kaiser,3},'SampleRate',Fs);

filterを使用してデータをフィルター処理します。フィルター出力の最初の 100 ポイントを、入力信号と同じ振幅と初期位相をもつ正弦波と重ね合わせてプロットします。

yout = filter(Hd,x); xin = cos(2*pi*100*t); plot(t,xin,t,yout) xlim([0 0.1]) xlabel('Time (s)') ylabel('Amplitude') legend('Input Sine Wave','Filtered Data') grid

フィルター処理されたデータの最初の 0.01 秒に着目すると、出力が入力に対して遅延していることがわかります。遅延は約 0.01 秒で、これはサンプルにおける FIR フィルターの長さのほぼ半分 $$(10\times 0.001)$$に相当します。

この遅延は、フィルターの位相応答によるものです。これらの例における FIR フィルターは、I 型の線形位相フィルターです。フィルターの群遅延は 10 サンプルです。

fvtoolを使用して群遅延をプロットします。

fvtool(Hd,'Analysis','grpdelay')

多くのアプリケーションで、位相の歪みは許容されています。このことは、位相応答が線形の場合に特に当てはまります。他のアプリケーションでは、ゼロ位相応答のフィルターの使用が望まれます。ゼロ位相応答は、非因果性フィルターでは技術的に不可能です。ただし、filtfiltによる因果性フィルターを使用してゼロ位相フィルター処理を実装することができます。

filtfiltを使用して入力信号をフィルター処理します。応答をプロットして、filterおよびfiltfiltで取得したフィルター出力を比較します。

yzp = filtfilt(Hd,x); plot(t,xin,t,yout,t,yzp) xlim([0 0.1]) xlabel('Time (s)') ylabel('Amplitude') legend('100-Hz Sine Wave','Filtered Signal','Zero-phase Filtering',...'Location','NorthEast')

上の図で、filtfilt出力には FIR フィルターの位相応答による遅延が現れないことがわかります。