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

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

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

定遅延フィルターの補正

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

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

Fs = 500;%采样率，单位为HzN = 500;%信号采样个数rng默认的；x =心电图(N) + 0.25 * randn (N, 1);%的波形t = (0: n - 1) / Fs;%的时间向量

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

Fnorm = 75 / (Fs / 2);%归一化频率df = designfilt (“lowpassfir”，“FilterOrder”, 70,“CutoffFrequency”, Fnorm);

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

grpdelay (df、2048 Fs)%标绘组延迟

D =意味着(grpdelay (df))采样滤波器延迟%

D = 35

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

y =过滤器(df, [x;0 (D, 1)]);%给输入数据追加D个零y = y (D + 1:结束);移动数据以补偿延迟情节(t t, x,, y,“r”，“线宽”, 1.5)标题(过滤后的波形的)包含(“时间(s)”)传说(“原始噪声信号”，“过滤信号”网格)在轴紧

周波数依存の遅延の補正

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

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

Fnorm = 75 / (Fs / 2);%归一化频率df = designfilt (“lowpassiir”，…“PassbandFrequency”Fnorm,…“FilterOrder”7…“PassbandRipple”, 1…“StopbandAttenuation”、60);

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

grpdelay (df, 2048,“一半”Fs)

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

日元=过滤器(df, x);%非线性相位滤波器-无延迟补偿y2 = filtfilt (df, x);%零相位执行-延迟补偿情节(t, x)在情节(t, y1,“r”，“线宽”1.5)情节(t, y2,“线宽”, 1.5)标题(过滤后的波形的)包含(“时间(s)”)传说(原始信号的，非线性相位IIR输出，…“零相位IIR输出”) xlim([0.25 0.55])网格在

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

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

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

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

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

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

Fs = 44100;y = audioread (“noisymusic.wav”）;

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

[P F] = pwelch (y)的(8192 1),8192/2,8192 Fs,“权力”）;helperFilterIntroductionPlot1 (F P [60 60], [-9.365 - -9.365],…｛“原始信号功率谱”，“60 Hz的语气”}）

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

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

《外交政策》= 1 e3;%通带频率，单位为Hz置= 1.4 e3;%阻带频率，单位为Hz美联社= 1;通频带纹波百分比，单位为dBAst = 95;%阻带衰减，单位为dBdf = designfilt (“lowpassfir”，“PassbandFrequency”《外交政策》,…“StopbandFrequency”置,“PassbandRipple”据美联社,,…“StopbandAttenuation”Ast,“SampleRate”Fs);

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

fvtool (df,“Fs”Fs,“FrequencyScale”，“日志”，…“FrequencyRange”，“指定freq.向量”，“FrequencyVector”F)

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

D =意味着(grpdelay (df));%过滤器延迟ylp =过滤器(df, [y;0 (D, 1)]);ylp = ylp (D + 1:结束);

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

Flp (Plp) = pwelch (ylp (8192 1), 8192/2, 8192 Fs,“权力”）;Flp helperFilterIntroductionPlot1 (F P、Plp…｛原始信号的，“低通滤过的信号”}）

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

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

Fs = f / 10;码= downsample (ylp 10);(Pds, Fds) = pwelch(码,(8192 1),8192/2,8192 Fs,“权力”）;helperFilterIntroductionPlot1 (F P Fds, Pds,…｛'信号以44100赫兹采样'，'下采样信号，Fs = 4410 Hz'}）

次に,IIRバンドストップフィルターで60 Hzトーンを除去します。阻止帯域が中心60 Hz,幅4赫兹になるようにします。シャープな周波数ノッチおよび小さい通過帯域リップルを比較的低い次数で実現できる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”，“日志”，…“FrequencyRange”，“指定freq.向量”，“FrequencyVector”Fds (Fds > F (2)))

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

yb = filtfilt (df、码);

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

yf =插值函数(yb, 10);Fs = f * 10;

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

[Pfinal, Ffinal] = pwelch (yf, (8192 1), 8192/2, 8192 Fs,“权力”）;helperFilterIntroductionPlot1 (F P Ffinal Pfinal,…｛原始信号的，“最终过滤信号”}）

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

要播放原始信号，取消下面两行注释% hplayer = audioplayer(y, Fs);% (hplayer)玩要播放降噪信号，取消下面两行注释% hplayer = audioplayer(yf, Fs);% (hplayer);

信号の微分

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

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

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

负载quakedrift.matFs = 1000;%采样率dt = 1 / f;%的时间差异t =(0:长度(漂移)1)* dt;%的时间向量

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

df = designfilt (“differentiatorfir”，“FilterOrder”, 50岁,…“PassbandFrequency”, 100,“StopbandFrequency”, 120,…“SampleRate”Fs);

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

HFVT = fvtool(df，[1 -1]，1, 1)“MagnitudeDisplay”，“零”，“Fs”Fs);传奇(hfvt50阶FIR差分器，“差函数的响应”）;

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

v1 = diff / dt(漂移);a1 = diff / dt (v1);v1 = [0;v1);a1 = [0;0;a1];

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

D =意味着(grpdelay (df));%过滤器延迟v2 =过滤器(df,漂移;0 (D, 1)]);v2 = v2 (D + 1:结束);a2 =过滤器(df, v2;0 (D, 1)]);a2 = a2 (D + 1:结束);v2 = v2 / dt;a2 = a2 / dt ^ 2;

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

helperFilterIntroductionPlot2 (t,漂移,v1、v2, a1, a2)

信号の積分

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

漏洩積分器を $$\mathrm{一个}＝0．999$$で使用します。漏洩積分器フィルターの振幅応答をプロットします。フィルターは,効果的に高周波ノイズを除去するローパスフィルターとして動作します。

fvtool (1 [1 -.999]“Fs”Fs)

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

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,漂移,d_leakyint、v_original v_leakyint)

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

まとめ

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

フィルター適用の詳細については,信号处理工具箱™ドキュメンテーションを参照してください。デジタルフィルターの設計方法の詳細については,実践に即したデジタルフィルター設計の紹介の例を参照してください。

参考文献

Proakis, J. G.和D. G. Manolakis。数字信号处理:原理、算法和应用englewood Cliffs，新泽西:Prentice-Hall出版社，1996年。

Orfanidis, s . J。信号处理概论englewood Cliffs，新泽西:Prentice-Hall出版社，1996年。

付録

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