Signal Processing Toolbox を使用したデータのフィルター処理

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

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

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

rng default

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

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

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

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

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

[h,f] = freqz(b,1,[],Fs);
plot(f,mag2db(abs(h)))
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)')
grid

信号にフィルターを適用し、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 フィルターの設計と実装を行う方法を示します。

コマンドラインで filterDesigner と入力してアプリを起動します。
[応答タイプ] を [ローパス] に設定します。
[設計法] を [FIR] に設定し、[ウィンドウ] 法を選択します。
[フィルター次数] で [次数指定] を選択します。次数を 20 に設定します。
[周波数仕様] で、[単位] を [Hz] に、[Fs] を 1000 に、[Fc] を 150 にそれぞれ設定します。

[フィルター設計] をクリックします。
[ファイル]、[エクスポート...] を選択して、FIR フィルターを係数またはフィルターオブジェクトとして MATLAB® ワークスペースにエクスポートします。この例では、フィルターをオブジェクトとしてエクスポートします。変数名を Hd と指定します。

[エクスポート] をクリックします。
コマンドウィンドウで、エクスポートしたフィルターオブジェクトを使用して入力信号をフィルター処理します。100 Hz の正弦波の最初の 10 周期の結果をプロットします。

y2 = filter(Hd,x);

plot(t,x,t,y2)
xlim([0 0.1])

xlabel('Time (s)')
ylabel('Amplitude')
legend('Original Signal','Filtered Data')

[ファイル]、[MATLAB コードを生成]、[フィルター設計関数] を選択して MATLAB 関数を生成し、自分の仕様でフィルターオブジェクトを作成します。

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

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

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

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

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

rng default

Fs = 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('Cycles/day')
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))
hold on
plot(w,unwrap(phiiir))
hold off

xlabel('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('Cycles/day')
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 を使って、信号の生成とローパスフィルターの設計を再度行います。これらの変数が既にワークプレースにある場合は、次のコードを実行する必要はありません。

rng default

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

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

% Using designfilt
Hd = 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 サンプルです。

群遅延をプロットします。

[gd,f] = grpdelay(Hd,[],Fs);

plot(f,gd)
xlabel('Frequency (Hz)')
ylabel('Group Delay (samples)')
grid

多くのアプリケーションで、位相の歪みは許容されています。このことは、位相応答が線形の場合に特に当てはまります。他のアプリケーションでは、ゼロ位相応答のフィルターの使用が望まれます。ゼロ位相応答は、非因果性フィルターでは技術的に不可能です。ただし、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 フィルターの位相応答による遅延が現れないことがわかります。

Signal Processing Toolbox を使用したデータのフィルター処理

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

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

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

ゼロ位相フィルター処理

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

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