Main Content

このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

実践に即したデジタル フィルター処理の紹介

この例では、デジタル フィルターを設計して、解析し、データのフィルター処理を行う方法を示します。これは次のような問題の解決に役立ちます。フィルターによって生じる遅延はどのように補正するのか、信号の歪みはどのようにして防ぐことができるのか、信号から望ましくない成分を除去するにはどういう方法があるのか、信号はどのようにして微分するのか、また積分するのか。

フィルターを使用することで、信号スペクトルを望ましい形にしたり、微分や積分のような数学演算を行ったりできます。以下において、フィルターが必要になったときにその利用が楽になるような、実践に即した概念をいくつか説明します。

この例ではデジタル フィルターの設計でなく、その適用に焦点を合わせて説明します。デジタル フィルターの設計方法に関する詳細については、実践に即したデジタル フィルター設計の紹介の例の例を参照してください。

フィルター処理によって生じる遅延の補正

デジタル フィルターを使用すると信号に遅延が生じます。フィルター特性によって、遅延が全周波数にわたって一定の場合と周波数に応じて変化する場合があります。遅延を補正するために行う処理は、遅延のタイプによって異なります。関数 grpdelay を使用すると、フィルター遅延を周波数関数としてとらえることができます。この関数出力を見ることによって、フィルター遅延が一定であるか、周波数に応じて変化する (つまり、遅延が周波数に依存する) かを特定することができます。

フィルター遅延が全周波数にわたって一定の場合、信号を時間シフトすることによって容易に補正できます。FIR フィルターでは通常、遅延が一定です。一方、遅延が周波数に応じて変化する場合、位相歪みが生じて信号波形が著しく変わることがあります。周波数に依存する遅延の補正は、遅延が一定の場合ほど簡単ではありません。IIR フィルターでは、周波数に依存する遅延が生じます。

定遅延フィルターの補正

前述のように、フィルターの群遅延を測定して周波数関数が一定であることを確認できます。関数 grpdelay を使用してフィルター遅延 D を測定し、入力信号に D 個の 0 を追加して出力信号を D サンプルだけ時間シフトすると、この遅延を補正できます。

ノイズを含んだ心電図信号について、75 Hz を上回る高周波ノイズを除去するためにフィルター処理する例を考えてみます。FIR ローパス フィルターをかけてから、フィルター遅延を補正します。これにより、ノイズを含んだ信号とフィルター処理後の信号が正しく整列され、比較できるようにプロットが互いに重ね合わせられます。

Fs = 500;                    % Sample rate in Hz
N = 500;                     % Number of signal samples
rng default;
x = ecg(N)'+0.25*randn(N,1); % Noisy waveform
t = (0:N-1)/Fs;              % Time vector

75 Hz のカットオフ周波数をもつ 70 次ローパス FIR フィルターを設計します。

Fnorm = 75/(Fs/2);           % Normalized frequency
df = designfilt('lowpassfir','FilterOrder',70,'CutoffFrequency',Fnorm);

フィルターの群遅延をプロットして、全周波数で群遅延が一定であり、このフィルターが線形位相であることを確認します。群遅延を使用してフィルター遅延を測定します。

grpdelay(df,2048,Fs)         % Plot group delay

Figure Figure 1: Group delay contains an axes object. The axes object with title Group delay, xlabel Frequency (Hz), ylabel Group delay (in samples) contains an object of type line.

D = mean(grpdelay(df)) % Filter delay in samples
D = 35

フィルター処理の前に、入力データ ベクトル x の最後に D 個の 0 を追加します。これにより、有用なサンプルはすべてフィルターから除外され、また、入力信号と遅延補正された出力信号が同じ長さになることが確実になります。データのフィルター処理を行い、出力信号を D サンプルだけ時間シフトして遅延を補正します。この最後の手順で、フィルターの過渡特性が効果的に除去されます。

y = filter(df,[x; zeros(D,1)]);  % Append D zeros to the input data
y = y(D+1:end);                  % Shift data to compensate for delay

plot(t,x,t,y,'linewidth',1.5)
title('Filtered Waveforms')
xlabel('Time (s)')
legend('Original Noisy Signal','Filtered Signal')
grid on
axis tight

Figure contains an axes object. The axes object with title Filtered Waveforms, xlabel Time (s) contains 2 objects of type line. These objects represent Original Noisy Signal, Filtered Signal.

周波数依存の遅延の補正

周波数に依存する遅延は信号の位相歪みを発生させます。このタイプの遅延を補正することは、遅延が一定である場合ほど簡単ではありません。オフライン処理が可能なアプリケーションの場合、関数 filtfilt を使用したゼロ位相フィルター処理によって周波数依存の遅延を除去できます。filtfilt は入力データの処理を順方向と逆方向の両方向に行うことにより、ゼロ位相フィルター処理を行います。主な効果はゼロ位相の歪みを得ること、つまり、データは 0 サンプルの定遅延をもつ同等のフィルターによってフィルター処理されます。その他の効果は、元のフィルターの伝達関数の振幅の 2 乗に等しいフィルター伝達関数と、元のフィルターの次数の 2 倍のフィルター次数が得られることです。

前節で定義した心電図信号を考えてみましょう。この信号を、遅延の補正がある場合とない場合に分けてフィルター処理します。75 Hz のカットオフ周波数をもつ 7 次ローパス IIR 楕円フィルターを設計します。

Fnorm = 75/(Fs/2); % Normalized frequency
df = designfilt('lowpassiir',...
               'PassbandFrequency',Fnorm,...
               'FilterOrder',7,...
               'PassbandRipple',1,...
               'StopbandAttenuation',60);

フィルターの群遅延をプロットします。群遅延が周波数に応じて変化し、フィルター遅延が周波数に依存しています。

grpdelay(df,2048,'half',Fs)

Figure Figure 2: Group delay contains an axes object. The axes object with title Group delay, xlabel Frequency (Hz), ylabel Group delay (in samples) contains an object of type line.

データをフィルター処理して、時間信号に対するそれぞれのフィルターの実装の効果を見てみます。ゼロ位相フィルター処理がフィルター遅延を効果的に除去しています。

y1 = filter(df,x);    % Nonlinear phase filter - no delay compensation
y2 = filtfilt(df,x);  % Zero-phase implementation - delay compensation

plot(t,x)
hold on
plot(t,y1,'r','linewidth',1.5)
plot(t,y2,'linewidth',1.5)
title('Filtered Waveforms')
xlabel('Time (s)')
legend('Original Signal','Non-linear phase IIR output',...
  'Zero-phase IIR output')
xlim([0.25 0.55])
grid on

Figure contains an axes object. The axes object with title Filtered Waveforms, xlabel Time (s) contains 3 objects of type line. These objects represent Original Signal, Non-linear phase IIR output, Zero-phase IIR output.

因果性のない順方向/逆方向のフィルター処理が可能で、フィルター応答を元の応答の 2 乗に変更できるアプリケーションの場合、ゼロ位相フィルター処理は優れたツールです。

定遅延が発生するフィルターは線形位相フィルターです。周波数依存の遅延が発生するフィルターは非線形位相フィルターです。

信号の望ましくないスペクトル成分の除去

フィルターは、信号からの望ましくないスペクトル成分の除去によく使用されます。これを行うため、さまざまな種類のフィルターからの選択が可能です。高周波数成分を除去する場合はローパス フィルターを選び、低周波数成分を除去する場合はハイパス フィルターを選びます。また、中間周波数帯域をそのまま残して低周波数成分も高周波数成分も除去する場合は、バンドパス フィルターを選びます。与えられた帯域の周波数成分を除去する場合は、バンドストップ フィルターを選びます。

電源供給ラインからのハムおよびホワイト ノイズを含むオーディオ信号について考えてみます。電源供給ラインからのハムは 60 Hz トーンが原因となって引き起こされます。ホワイト ノイズはオーディオ帯域全体にわたって存在する信号です。

オーディオ信号を読み込みます。サンプル レートを 44.1 kHz に指定します。

Fs = 44100;
y = audioread('noisymusic.wav');

信号のパワー スペクトルをプロットします。赤い三角形のマーカーは、オーディオ信号に干渉する強い 60 Hz トーンを示します。

[P,F] = pwelch(y,ones(8192,1),8192/2,8192,Fs,'power');
helperFilterIntroductionPlot1(F,P,[60 60],[-9.365 -9.365],...
  {'Original signal power spectrum', '60 Hz Tone'})

Figure contains an axes object. The axes object with xlabel Frequency in kHz, ylabel Power Spectrum (dB) contains 2 objects of type line. One or more of the lines displays its values using only markers These objects represent Original signal power spectrum, 60 Hz Tone.

最初に、できるだけ多くのホワイト ノイズ スペクトル成分をローパス フィルターで除去します。フィルターの通過帯域の値は、ノイズ除去と高周波数成分の損失によるオーディオ品質の劣化との適切なトレードオフが得られる設定にすべきです。60 Hz ハムを除去する前にローパス フィルターをかけると、帯域制限された信号のダウンサンプリングが可能になり、非常に便利です。信号のレートを低くすると、よりシャープで狭い 60 Hz バンドストップ フィルターを低いフィルター次数で設計できます。

通過帯域周波数 1 kHz、阻止帯域周波数 1.4 kHz のローパス フィルターを設計します。最小次数設計を選択します。

Fp = 1e3;    % Passband frequency in Hz
Fst = 1.4e3; % Stopband frequency in Hz
Ap = 1;      % Passband ripple in dB
Ast = 95;    % Stopband attenuation in dB

df = designfilt('lowpassfir','PassbandFrequency',Fp,...
                'StopbandFrequency',Fst,'PassbandRipple',Ap,...
                'StopbandAttenuation',Ast,'SampleRate',Fs);

フィルター応答を解析します。

fvtool(df,'Fs',Fs,'FrequencyScale','log',...
  'FrequencyRange','Specify freq. vector','FrequencyVector',F)

Figure Figure 3: Magnitude Response (dB) contains an axes object. The axes object with title Magnitude Response (dB), xlabel Frequency (kHz), ylabel Magnitude (dB) contains 2 objects of type line.

データをフィルター処理し、遅延を補正します。

D = mean(grpdelay(df)); % Filter delay
ylp = filter(df,[y; zeros(D,1)]);
ylp = ylp(D+1:end);

ローパス フィルター処理を行った信号のスペクトルを見てみます。1400 Hz を超える周波数成分が除去されています。

[Plp,Flp] = pwelch(ylp,ones(8192,1),8192/2,8192,Fs,'power');
helperFilterIntroductionPlot1(F,P,Flp,Plp,...
  {'Original signal','Lowpass filtered signal'})

Figure contains an axes object. The axes object with xlabel Frequency in kHz, ylabel Power Spectrum (dB) contains 2 objects of type line. These objects represent Original signal, Lowpass filtered signal.

上のパワー スペクトルのプロットより、ローパス フィルター処理を行った信号の無視できない周波数成分の最大値は 1400 Hz であることがわかります。サンプリング理論により、サンプル レート 2×1400=2800 は信号を正確に表すのに十分です。しかし、必要以上のサンプルを処理する必要があって無駄な 44100 Hz のサンプル レートを使用しています。信号をダウンサンプリングしてサンプル レートを下げ、処理サンプル数を減らすことで計算負荷を軽減できます。サンプル レートを低くすると、60 Hz ノイズ除去に必要なよりシャープで狭いバンドストップ フィルターを低いフィルター次数で設計することができます。

サンプル レート Fs/10 = 4.41 kHz を得るために、ローパス フィルター処理を行った信号を係数 10 でダウンサンプリングします。ダウンサンプリング前後の信号のスペクトルをプロットします。

Fs = Fs/10;
yds = downsample(ylp,10);

[Pds,Fds] = pwelch(yds,ones(8192,1),8192/2,8192,Fs,'power');
helperFilterIntroductionPlot1(F,P,Fds,Pds,...
  {'Signal sampled at 44100 Hz', 'Downsampled signal, Fs = 4410 Hz'})

Figure contains an axes object. The axes object with xlabel Frequency in kHz, ylabel Power Spectrum (dB) contains 2 objects of type line. These objects represent Signal sampled at 44100 Hz, Downsampled signal, Fs = 4410 Hz.

次に、IIR バンドストップ フィルターで 60 Hz トーンを除去します。阻止帯域が中心 60 Hz、幅 4 Hz になるようにします。シャープな周波数ノッチおよび小さい通過帯域リップルを比較的低い次数で実現できる IIR フィルターを選択します。

df = designfilt('bandstopiir','PassbandFrequency1',55,...
               'StopbandFrequency1',58,'StopbandFrequency2',62,...
               'PassbandFrequency2',65,'PassbandRipple1',1,...
               'StopbandAttenuation',60,'PassbandRipple2',1,...
               'SampleRate',Fs,'DesignMethod','ellip');                          

振幅応答を解析します。

fvtool(df,'Fs',Fs,'FrequencyScale','log',...
  'FrequencyRange','Specify freq. vector','FrequencyVector',Fds(Fds>F(2)))

Figure Figure 4: Magnitude Response (dB) contains an axes object. The axes object with title Magnitude Response (dB), xlabel Frequency (kHz), ylabel Magnitude (dB) contains 2 objects of type line.

位相歪みを防ぐために、ゼロ位相フィルター処理を実行します。関数 filtfilt を使用してデータを処理します。

ybs = filtfilt(df,yds);

最後に、信号をアップサンプリングして、オーディオ サウンド カードに準拠した元のオーディオ サンプル レート 44.1 kHz に戻します。

yf = interp(ybs,10);
Fs = Fs*10;

最後に、元の信号のスペクトルとフィルター処理後のスペクトルを見てみます。フィルターによって、高周波ノイズ フロアと 60 Hz トーンが減衰しています。

[Pfinal,Ffinal] = pwelch(yf,ones(8192,1),8192/2,8192,Fs,'power');
helperFilterIntroductionPlot1(F,P,Ffinal,Pfinal,...
  {'Original signal','Final filtered signal'})

Figure contains an axes object. The axes object with xlabel Frequency in kHz, ylabel Power Spectrum (dB) contains 2 objects of type line. These objects represent Original signal, Final filtered signal.

コンピューターにサウンドカードが搭載されているなら、処理前後の信号を再生できます。上述したように、最終的にオーディオ ファイルの 60 Hz ハムと高周波ノイズを効果的に減衰させることができました。

% To play the original signal, uncomment the next two lines
% hplayer = audioplayer(y, Fs);
% play(hplayer)

% To play the noise-reduced signal, uncomment the next two lines
% hplayer = audioplayer(yf, Fs);
% play(hplayer);

信号の微分

MATLAB 関数 diff による信号の微分には、出力時ノイズ レベル増大の可能性という欠点があります。このため、より良い選択として微分器フィルターを使用します。微分器フィルターは、目的の帯域で微分器として、他の全周波数では減衰器として動作し、高周波ノイズを効果的に除去します。

例として、地震発生時の建物の床の変位速度を解析します。変位またはドリフトの測定値は、地震の条件下で 3 階建ての試験構造物の 1 階で記録されて quakedrift.mat ファイルに保存されました。データ ベクトルの長さは 10e3、サンプル レートは 1 kHz、測定単位は cm です。

変位データを微分して、地震発生時の建物の床の速度と加速度の推定値を求めます。diff および FIR 微分器フィルターを使用した結果を比較します。

load quakedrift.mat 

Fs  = 1000;                 % Sample rate
dt = 1/Fs;                  % Time differential
t = (0:length(drift)-1)*dt; % Time vector

ほとんどの信号エネルギーが存在する帯域である 100 Hz の通過帯域周波数の 50 次微分器フィルターを設計します。フィルターの阻止帯域周波数を 120 Hz に設定します。

df = designfilt('differentiatorfir','FilterOrder',50,...
                'PassbandFrequency',100,'StopbandFrequency',120,...
                'SampleRate',Fs);

関数 diff は、応答 H(Z)=1-Z-1 をもつ 1 次 FIR フィルターと見なすことができます。50 次微分器 FIR フィルターの振幅応答と関数 diff の応答を比較するために、FVTool を使用します。明らかに、0 ~ 100 Hz の通過帯域領域で両方の応答が等価です。しかし、阻止帯域領域では、50 次微分器フィルターが成分を減衰させているのに対し、diff の応答は成分を増幅させています。これが実質的に、高周波ノイズのレベルを上げています。

hfvt = fvtool(df,[1 -1],1,'MagnitudeDisplay','zero-phase','Fs',Fs);
legend(hfvt,'50th order FIR differentiator','Response of diff function');

Figure Figure 5: Zero-phase Response contains an axes object. The axes object with title Zero-phase Response, xlabel Frequency (Hz), ylabel Amplitude contains 2 objects of type line. These objects represent 50th order FIR differentiator, Response of diff function.

関数 diff を使用して微分します。diff の演算によって欠損しているサンプルを補正するために、0 を追加します。

v1 = diff(drift)/dt;
a1 = diff(v1)/dt;

v1 = [0; v1];
a1 = [0; 0; a1];

50 次微分器 FIR フィルターを使用して微分し、遅延を補正します。

D = mean(grpdelay(df)); % Filter delay
v2 = filter(df,[drift; zeros(D,1)]);
v2 = v2(D+1:end);
a2 = filter(df,[v2; zeros(D,1)]);
a2 = a2(D+1:end);
v2 = v2/dt;
a2 = a2/dt^2;

床の変位データ点をいくつかプロットします。diff および 50 次微分器 FIR フィルターによって算出した速度と加速度のデータ点もいくつかプロットします。diff を使用した場合、速度の推定値ではノイズがわずかに増幅し、加速度の推定値ではノイズが大幅に増幅していることがわかります。

helperFilterIntroductionPlot2(t,drift,v1,v2,a1,a2)

Figure contains an axes object. The axes object with xlabel Time in seconds, ylabel Displacement (cm) contains an object of type line.

Figure contains an axes object. The axes object with xlabel Time in seconds, ylabel speed (cm/s) contains 2 objects of type line. These objects represent Estimated speed using diff, Estimated speed using FIR filter.

Figure contains an axes object. The axes object with xlabel Time in seconds, ylabel acceleration (cm/s toThePowerOf 2) baseline contains 2 objects of type line. These objects represent Estimated acceleration using diff, Estimated acceleration using FIR filter.

信号の積分

漏洩積分器フィルターは、伝達関数 H(Z)=1/[1-cZ-1] を使用する全極フィルターです。c は定数で、フィルターの安定性を確保するためには 1 未満でなければなりません。当然ですが、c が 1 に近づくにつれて、漏洩積分器が diff の伝達関数の逆関数に近づきます。前節で求めた加速度と速度の推定値を漏洩積分器で積分して、速度とドリフトをそれぞれ取り戻します。関数 diff を使用して求めた推定値のほうがノイズが多いので、こちらを使用します。

漏洩積分器を a=0.999 で使用します。漏洩積分器フィルターの振幅応答をプロットします。フィルターは、効果的に高周波ノイズを除去するローパス フィルターとして動作します。

fvtool(1,[1 -.999],'Fs',Fs)

Figure Figure 6: Magnitude Response (dB) contains an axes object. The axes object with title Magnitude Response (dB), xlabel Frequency (Hz), ylabel Magnitude (dB) contains an object of type line.

漏洩積分器で、速度と加速度をフィルター処理します。時間差で乗算します。

v_original = v1;
a_original = a1;

d_leakyint = filter(1,[1 -0.999],v_original);
v_leakyint = filter(1,[1 -0.999],a_original);

d_leakyint = d_leakyint * dt;
v_leakyint = v_leakyint * dt;

変位と速度の推定値をプロットし、元の信号と比較します。

helperFilterIntroductionPlot3(t,drift,d_leakyint,v_original,v_leakyint)

Figure contains 2 axes objects. Axes object 1 with ylabel Displacement (cm) contains 2 objects of type line. These objects represent Leaky integrator estimate, Original displacement. Axes object 2 with xlabel Time in seconds, ylabel Speed (cm/s) contains 2 objects of type line. These objects represent Leaky integrator estimate, Original speed.

関数 cumsum および cumtrapz を使用しても信号を積分できます。結果は、漏洩積分器で求めた結果と同様になります。

まとめ

この例では、線形位相フィルターおよび非線形位相フィルターについて、また、それぞれのフィルター タイプによって生じる位相遅延の補正方法について説明しました。さらに、フィルターをかけて望ましくない周波数成分を信号から取り除く方法やローパス フィルターで帯域を制限した後で信号をダウンサンプリングする方法について説明しました。最後に、デジタル フィルター設計を使用した、信号の微分および積分方法を説明しました。この例の全体を通して、解析ツールを使用してフィルターの応答と群遅延を確認する方法についても説明しました。

フィルター適用の詳細については、Signal Processing Toolbox™ ドキュメンテーションを参照してください。デジタル フィルターの設計方法の詳細については、実践に即したデジタル フィルター設計の紹介の例を参照してください。

参考文献

Proakis, J. G., and D. G. Manolakis.Digital Signal Processing: Principles, Algorithms, and Applications.Englewood Cliffs, NJ: Prentice-Hall, 1996.

Orfanidis, S. J. Introduction To Signal Processing.Englewood Cliffs, NJ: Prentice-Hall, 1996.

付録

この例では、次の補助関数が使用されています。

参考

| | | | | | | |