数字滤波实用导论

打开生活的脚本

这个例子展示了如何设计、分析和应用数字滤波器到您的数据。它将帮助您回答以下问题:我如何补偿滤波器引入的延迟?、如何避免信号失真?，How do I remove unwanted content from my signal?, How do I differentiate my signal?, and How do I integrate my signal?

滤波器可用于以所需方式塑造信号频谱，或执行微分和积分等数学运算。在下面的内容中，您将学习一些实用的概念，这些概念将在您需要时简化过滤器的使用。

这个例子着重于数字滤波器的应用而不是其设计。如果你想了解更多关于如何设计数字滤波器，请参阅数字滤波器设计实用介绍实例

滤波引起的延迟补偿

数字滤波器在信号中引入延迟。根据滤波器的特性，延迟可以在所有频率上是恒定的，也可以随频率而变化。延迟的类型决定了你必须采取什么行动来弥补它。的grpdelay函数允许您将滤波器延迟视为频率的函数。通过查看此函数的输出，可以确定滤波器的延迟是恒定的，还是随频率变化（换句话说，它是否与频率相关）。

滤波器延迟在所有频率上都是恒定的，可以很容易地通过在时间上移动信号来补偿。FIR滤波器通常有恒定的延迟。另一方面，随频率变化的延迟会引起相位失真，并显著改变信号波形。对频率相关延迟的补偿不像对常量延迟的补偿那么简单。IIR滤波器引入频率相关延迟。

对恒定滤波器延迟的补偿

如前所述，可以测量滤波器的延迟组，以验证它是频率的常数函数。你可以使用grpdelay函数用于测量滤波器延迟D，并通过向输入信号添加D零并将输出信号按时间偏移D个采样来补偿该延迟。

考虑一个有噪声的心电图信号，您想要过滤掉高于75hz的高频噪声。你想要应用一个FIR低通滤波器和补偿滤波器延迟，以便噪声和滤波信号正确对齐，并可以绘制在彼此之上进行比较。

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

设计一个截止频率为75hz的70阶低通FIR滤波器。

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

绘制滤波器的组延迟，以验证它是常数跨越所有频率，表明滤波器是线性相位。使用组延迟来测量滤波器的延迟。

grpdelay (df、2048 Fs)图群延迟

Figure Filter Visualization Tool-Group delay包含一个轴对象和uitoolbar、uimenu类型的其他对象。标题为Group delay的轴对象包含一个line类型的对象。

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

D = 35

在滤波之前，在输入数据向量x的末尾添加D个零。这确保了所有有用的样本都从滤波器中清除出来，并且输入信号和延迟补偿输出信号具有相同的长度。对数据进行滤波，并通过D采样对输出信号进行移位来补偿延时。最后一步有效地删除了瞬态滤波器。

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

图中包含一个轴对象。标题为“过滤波形”的轴对象包含两个类型为line的对象。这些对象分别代表原始噪声信号和滤波信号。

补偿频率相关延迟

频率相关的延迟导致信号的相位失真。这种类型的延迟的补偿不像常量延迟的补偿那么简单。如果您的应用程序允许脱机处理，您可以通过使用filtfilt函数。filtfilt通过对输入数据进行正向和反向处理，实现零相位滤波。主要的效果是获得零相位失真，也就是说，你用一个等效的滤波器过滤数据，它有一个恒定的0采样延迟。其他的效果是你得到一个滤波器传递函数它等于原滤波器传递函数的平方大小，一个滤波器的阶数是原滤波器阶数的两倍。

考虑上一节定义的心电信号。滤波器这个信号与没有延迟补偿。设计一个截止频率为75hz的七阶低通IIR椭圆滤波器。

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

绘制滤波器的组延迟。组延迟随频率变化，说明滤波器延迟与频率有关。

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

Figure Filter Visualization Tool-Group delay包含一个轴对象和uitoolbar、uimenu类型的其他对象。标题为Group delay的轴对象包含一个line类型的对象。

过滤数据，并查看每个滤波器实现对时间信号的影响。零相位滤波有效地消除了滤波器延迟。

y1=滤波器（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])网格在

图中包含一个轴对象。标题为“过滤波形”的轴对象包含3个类型为line的对象。这些对象分别代表原始信号、非线性相位IIR输出、零相位IIR输出。

如果您的应用程序允许非因果的前向/后向过滤操作，并且允许将过滤器响应更改为原始响应的平方，那么零相位过滤是一个很好的工具。

引入恒定延迟的滤波器是线性相位滤波器。引入频率相关延迟的滤波器是非线性相位滤波器。

去除信号中不需要的光谱内容

滤波器通常用于去除信号中不需要的光谱内容。您可以从各种过滤器中进行选择。当你想要去除高频内容时，你可以选择低通滤波器，或者当你想要去除低频内容时，你可以选择高通滤波器。您也可以选择带通滤波器，以去除低频和高频内容，而留下一个中间频带的频率完好无损。当你想要去除给定频带上的频率时，你可以选择带阻滤波器。

考虑具有电源线嗡嗡声和白噪声的音频信号。电源线嗡嗡声是由60 Hz的音调引起的。白噪声是一种存在于所有音频带宽上的信号。

加载音频信号。指定44.1 kHz的采样率。

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

绘制信号的功率谱。红色三角形标记显示强烈的60hz音调干扰音频信号。

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

图中包含一个轴对象。轴对象包含两个类型为line的对象。这些对象代表原始信号功率谱，60hz音调。

首先可以使用低通滤波器去除尽可能多的白噪声频谱内容。滤波器的通频带应该设置为一个值，在噪声降低和由于高频内容丢失而导致的音频退化之间提供一个良好的平衡。在去除60hz的嗡嗡声之前应用低通滤波器是非常方便的，因为你将能够向下采样带限信号。较低的速率信号将允许您设计一个更尖锐和更窄的60赫兹带阻滤波器较小的滤波器订单。

设计一个通频带频率为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,“频率标度”，“日志”，．..“频率范围”，“指定freq.向量”，“FrequencyVector”F)

Figure Filter Visualization Tool-震级响应（dB）包含一个轴对象和uitoolbar、uimenu类型的其他对象。标题为震级响应（dB）的轴对象包含两个line类型的对象。

过滤数据并补偿延迟。

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

查看低通滤波信号的频谱。1400 Hz以上的频率内容已被删除。

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

图中包含一个轴对象。轴对象包含两个类型为line的对象。这些对象代表原始信号，低通滤波信号。

从上面的功率谱图可以看出，低通滤波信号的最大不可忽略的频率含量是在1400 Hz。根据抽样定理，抽样率为 $2 \times 1400 ＝ 2800$ Hz足以正确表示信号，但是，您使用的采样率为44100 Hz，这是一种浪费，因为您将需要处理比所需更多的样本。您可以通过减少需要处理的采样数来降低信号采样率和计算负载。较低的采样率还将允许您设计更清晰、更窄的带阻滤波器，以较小的滤波器阶数消除60 Hz噪声。

对低通滤波信号进行10倍的下采样，以获得Fs/10 = 4.41 kHz的采样率。绘制下采样前后的信号频谱。

Fs = f / 10;码= downsample (ylp 10);(Pds, Fds) = pwelch(码,(8192 1),8192/2,8192 Fs,“权力”); helperFilterIntroductionPlot1（F、P、Fds、Pds、，．..｛“信号以44100hz采样”，'下采样信号，Fs = 4410hz '})

图中包含一个Axis对象。Axis对象包含2个line类型的对象。这些对象表示在44100 Hz下采样的信号，下采样信号，Fs=4410 Hz。

现在使用IIR带阻滤波器去除60hz的音调。让阻带宽度为4hz，以60hz为中心。选择一个IIR滤波器，以达到尖锐的频率陷波，小通带纹波和相对较低的阶数。

df = designfilt (“bandstopiir”，“通带频率1”现年55岁的．..“StopbandFrequency1”58岁的“StopbandFrequency2”, 62,．..“通带频率2”, 65,“PassbandRipple1”1.．..“StopbandAttenuation”现年60岁的“PassbandRipple2”1.．..“SampleRate”Fs,“设计方法”，“ellip”）;

分析震级响应。

fvtool（df，“财政司司长”Fs,“频率标度”，“日志”，．..“频率范围”，“指定freq.向量”，“FrequencyVector”Fds (Fds > F (2)))

Figure Filter Visualization Tool-震级响应（dB）包含一个轴对象和uitoolbar、uimenu类型的其他对象。标题为震级响应（dB）的轴对象包含两个line类型的对象。

执行零相位滤波以避免相位失真。使用filtfilt函数来处理数据。

yb = filtfilt (df、码);

最后，对信号进行上采样，使其恢复到与声卡兼容的原始音频采样率44.1 kHz。

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

最后看一下原始和处理过的信号的频谱。高频底噪声和60赫兹的音调已被滤波器衰减。

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

图中包含一个轴对象。轴对象包含两个类型为line的对象。这些对象代表原始信号，最终滤波信号。

如果你的电脑有声卡，你可以听信号处理前后。如上所述，最终结果是有效地衰减了音频文件中的60hz嗡嗡声和高频噪声。

%要播放原始信号，请取消对下两行的注释%hplayer=音频播放器（y，Fs）；% (hplayer)玩要播放降噪信号，取消注释下面两行% hplayer = audioplayer(yf, Fs);%播放（hplayer）；

区分一个信号

MATLABdiff函数区分信号的缺点是，您可能会在输出时增加噪声水平。一个更好的选择是使用微分器滤波器，在感兴趣的频带中作为微分器，在所有其他频率上作为衰减器，有效地去除高频噪声。

作为一个例子，分析地震期间建筑地板的位移速度。在地震条件下，在三层测试结构的第一层记录位移或漂移测量值，并保存在quakedrift.mat文件中。数据向量的长度为10e3，采样率为1 kHz，测量单位为厘米。

区分位移数据，以获得地震期间建筑物地板的速度和加速度估计值。使用diff和FIR微分器滤波器比较结果。

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

设计一个50阶微分器滤波器，其通频带频率为100hz，这是信号能量的主要来源。设置滤波器的阻带频率为120hz。

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

的diff函数可以看作是一个具有响应的一阶FIR滤波器 $H （ Z) ＝ 1 - Z^{- 1}$ ．使用FVTool比较50阶微分器FIR滤波器的幅值响应与该滤波器的响应diff函数。显然，这两个响应在通带区域(从0到100hz)是等效的。然而，在阻带区域，50阶滤波器衰减分量，而差分响应放大分量。这有效地增加了高频噪声的水平。

HFVT = fvtool(df，[1 -1]，1, 0);“震级显示”，“零相位”，“财政司司长”Fs);传奇(hfvt“50阶FIR微分器”，“差分函数的响应”）;

Figure Filter Visualization Tool-Zero phase Response包含一个轴对象和uitoolbar、uimenu类型的其他对象。标题为Zero phase Response的轴对象包含2个line类型的对象。这些对象表示50阶FIR微分器，即diff函数的响应。

区分使用diff函数。添加零来补偿由于diff操作而丢失的样本。

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

使用50阶FIR滤波器进行微分和延时补偿。

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)

图中包含一个轴对象。axis对象包含一个类型为line的对象。

图中包含一个axes对象。axes对象包含2个line类型的对象。这些对象表示使用diff的估计速度，使用FIR滤波器的估计速度。

图中包含一个轴对象。轴对象包含两个类型为line的对象。这些对象表示估计加速度使用diff，估计加速度使用FIR滤波器。

整合一个信号

漏式积分器滤波器是一种具有传递函数的全极滤波器 $H （ Z) ＝ 1 / ［ 1 - c Z^{- 1} ］$ 在哪里 $c$ 为必须小于1的常数，以确保滤波器的稳定性。这是不足为奇的 $c$ 趋近于1，漏积器趋近于diff传递函数。将漏积器应用于前一节得到的加速度和速度估计，分别得到速度和漂移。使用得到的估计的diff因为它们噪音更大。

使用漏泄积分器 $一个＝ 0 ． 999$ ．绘制漏出的积分器滤波器的幅度响应。该滤波器作为低通滤波器有效地消除高频噪声。

fvtool（1[1-.999]，“财政司司长”Fs)

Figure Filter Visualization Tool-震级响应（dB）包含轴对象和uitoolbar、uimenu类型的其他对象。标题为震级响应（dB）的轴对象包含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,漂移,d_leakyint、v_original v_leakyint)

图中包含2个轴对象。坐标轴对象1包含2个类型为line的对象。这些对象代表Leaky积分估计，原始位移。axis对象2包含2个类型为line的对象。这些对象表示漏积分估计，原始速度。

你也可以用cumsum和cumtrapz功能。结果将与使用泄漏积分器获得的结果相似。

结论

在这个例子中，您学习了线性和非线性相位滤波器，并学习了如何补偿由每种滤波器类型引入的相位延迟。您还学习了如何应用滤波器从信号中去除不需要的频率成分，以及如何在用低通滤波器限制其带宽后对信号进行下采样。最后，您学习了如何使用数字滤波器设计来区分和集成信号。在整个示例中，您还学习了如何使用分析工具查看过滤器的响应和组延迟。

有关过滤器应用程序的更多信息，请参阅信号处理工具箱™文档。有关如何设计数字滤波器的更多信息，请参阅数字滤波器设计实用介绍实例