数字滤波的实用介绍

打开脚本

这个示例展示了如何设计、分析和对数据应用数字过滤器。它将帮助您回答这样的问题:我如何补偿由过滤器引起的延迟?，如何避免信号失真?,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样品中时移的输出信号补偿该延迟。

考虑一个有噪声的心电图信号，您希望对其进行滤波以去除75赫兹以上的高频噪声。你想应用一个FIR低通滤波器，并补偿滤波器的延迟，使噪声和滤波后的信号正确地对齐，并可以绘制在彼此的上方进行比较。

FS = 500;以Hz％采样率N = 500;信号样本的数量％RNG默认的;x =心电图(N) + 0.25 * randn (N, 1);%的波形t = (0: n - 1) / Fs;％的时间矢量%设计70阶低通FIR滤波器，截止频率为75hz。Fnorm = 75 / (Fs / 2);%归一化频率df = designfilt (“lowpassfir”,“FilterOrder”,70,“CutoffFrequency”，FNORM）;

画出滤波器的群延迟，以确认它是横跨表明过滤器是线性相位的所有频率恒定。使用群时延测量滤波器的延迟。

grpdelay（DF，2048，FS）图群延迟d =平均值（grpdelay（DF））在样品％滤波器延迟

D = 35

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

y =过滤器(df, [x;0 (D, 1)]);％追加d零到输入数据Y = Y（d + 1：结束）;移位数据以补偿延迟图绘制(t t, x,, y,“r”,'行宽',1.5);标题(“滤波的波形”);包含('时间（s）')传说(“原始噪声信号”,“滤波信号”);网格上轴紧

补偿频率相关的延迟

频率相关的延迟会导致信号的相位失真。补偿这种类型的延迟不像补偿常数延迟那么简单。如果应用程序允许脱机处理，则可以通过使用filtfilt函数。filtfilt通过处理正向和反向的输入数据来执行零相位滤波。主要的影响是你获得零相位失真，即。,you filter data with an equivalent filter that has a constant delay of 0 samples. Other effects are that you get a filter transfer function which equals the squared magnitude of the original filter transfer function, and a filter order that is double the order of the original filter.

考虑前面部分定义的心电信号。滤波这个信号与没有延迟补偿。

％设计一个7阶低通IIR椭圆滤波器的截止频率750hz的%。Fnorm = 75 / (Fs / 2);%归一化频率df = designfilt ('lowpassiir',...'PassbandFrequency'，FNORM，...“FilterOrder”7...“PassbandRipple”1，...'StopbandAttenuation'，60）;

绘制滤波器的群时延，注意它随频率的变化，表明滤波器时延与频率有关。

grpdelay（DF，2048，“一半”，FS）

在时间信号的每个过滤器实现的效果滤镜中的数据和外观。

Y1 =滤波器（DF，X）;％的非线性相位滤波器 - 没有延迟补偿y2 = filtfilt (df, x);%零相位执行-延迟补偿图绘制(t, x);持有上情节(t, y1,“r”,'行宽',1.5);图（T，Y2，'行宽',1.5);标题(“滤波的波形”);包含('时间（s）')传说(“原始信号”,“非线性相位IIR输出”,...“零相位IIR输出”);ax =轴;轴([0.25 0.55 ax(3:4)])网格上

通知零相有效过滤如何移除滤波器延迟。

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

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

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

滤波器通常用于去除来自不想要的信号频谱内容。您可以从各种过滤器要做到这一点选择。当你想要的时候，你要删除的低频成分去除高频内容，或者高通滤波器您可以选择一个低通滤波器。也可以选择一个带通滤波器，以除去低和高频率成分，同时留下的频率的中间带完整无缺。当你想在给定的频段频率删除您选择一个带阻滤波器。

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

加载音频信号。

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赫兹音”}）

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

设计通带频率为1khz，阻带频率为1.4 kHz的低通滤波器。选择最小订单设计。

《外交政策》= 1 e3;％通带频率以Hz置= 1.4 e3;％阻带频率以HzAP = 1;%通带纹波在分贝AST = 95;％阻带衰减，以dB%设计过滤器df = designfilt (“lowpassfir”,'PassbandFrequency'《外交政策》,...“StopbandFrequency”置,“PassbandRipple”，鸭，...'StopbandAttenuation'Ast,“SampleRate”Fs);％分析滤波器响应hfvt = fvtool（DF，'FS'，FS，'FrequencyScale',“日志”,...'频率范围',“指定频率。向量','FrequencyVector'F);

％过滤数据和补偿延迟D =意味着(grpdelay (df));%过滤器延迟YLP =滤波器（DF，[Y;零（d，1）]）;YLP = YLP（d + 1：结束）;接近（hfvt）

看看低通滤波信号的频谱。请注意，频率内容超过1400赫兹已被删除。

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

从上面的功率谱图可以看出，低通滤波信号的最大不可忽略的频率含量是在1400 Hz。根据采样定理，2*1400 = 2800 Hz的采样频率就足以正确地表示信号，但是，您使用的是44100 Hz的采样率，这是一种浪费，因为您将需要处理更多的样本。您可以对信号进行降采样，以减少采样率，并通过减少需要处理的样本数量来减少计算负载。一个较低的采样率也将允许你设计一个更尖锐和更窄的带阻滤波器，需要消除60赫兹的噪音，一个更小的滤波器阶。

下采样由10倍的低通滤波后的信号以获得FS / 10 = 4.41千赫的采样率。之前和之后采样的绘制信号的频谱。

FS = FS / 10;YDS =下采样（YLP，10）;[PDS，FDS] = pwelch（YDS，酮（8192,1），8192 / 2,8192，FS“权力”);helperFilterIntroductionPlot1（F，P，FDS，PDS，...{“在44100hz采样的信号”,'下采样信号，FS = 4410赫兹'}）

现在，除去使用IIR带阻滤波器的60 Hz的音调。让阻带具有4赫兹，60赫兹为中心的宽度。我们选择IIR滤波器以获得清晰的频率缺口，小通带波纹，以及相对低的顺序。使用处理数据filtfilt避免相位失真。

%设计过滤器df = designfilt ('bandstopiir','PassbandFrequency1'55，...'StopbandFrequency1'58岁的'StopbandFrequency2'62，...'PassbandFrequency2'65，“PassbandRipple1”1，...'StopbandAttenuation'60，“PassbandRipple2”1，...“SampleRate”，FS，“DesignMethod”,“椭球”);％分析幅度响应hfvt = fvtool（DF，'FS'，FS，'FrequencyScale',“日志”,...'频率范围',“指定频率。向量','FrequencyVector'，FDS（FDS> F（2）））;

执行零相位滤波以避免失真。

yb = filtfilt (df、码);

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

YF = interp的（YBS，10）;FS = FS * 10;

最后看一下原始信号和处理过的信号的频谱。注意高频噪声地板和60hz的音调是如何被滤波器衰减的。

[Pfinal, Ffinal] = pwelch (yf, (8192 1), 8192/2, 8192 Fs,“权力”);关闭（hfvt）helperFilterIntroductionPlot1（F，P，Ffinal，Pfinal，...{原始信号的,“最终滤波信号”}）

听信号前后处理。如上所述，最终的结果是您已经有效地减弱了60赫兹的嗡嗡声和音频文件上的高频噪声。

％播放原始信号hplayer = audioplayer(y, Fs);玩(hplayer);播放降噪信号hplayer = audioplayer(yf, Fs);玩(hplayer);

区分信号

在MATLABdiff功能划分与缺点，即可以潜在地增加的噪声电平在输出的信号。一个更好的选择是使用微分滤波器充当在感兴趣的频带中的微分器，并且如在所有其他频率的衰减器，有效地去除高频噪声。

以某建筑物为例，分析其在地震作用下的位移速度。在地震条件下，位移或漂移测量记录在一个三层测试结构的一楼，并保存在地震裂缝中。垫文件。数据向量长度为10e3，采样率为1khz，测量单位为cm。

差异化的位移数据在地震中获得的速度和建筑物楼板的加速度的估计。使用比较差异和FIR滤波器区分的结果。

负载quakedrift.matFs = 1000;%采样率dt = 1 / f;％的时间微分T =（0：长度（漂移）-1）* dt的;％的时间矢量

设计一个通频带频率为100hz的50阶微分器滤波器，这是发现大部分信号能量的带宽。将滤波器的阻带频率设置为120hz。

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

的diff功能可以被看作是与响应的第一阶FIR滤波器 $$H(Z) = 1 - Z^{-1$ 。用FVTool比较了第50阶微分器FIR滤波器的幅值响应和第50阶微分器FIR滤波器的响应diff函数。显然，这两种反应是在通带区域（从0到100Hz）等效。然而，在阻带区域，第50位顺序滤波器衰减而DIFF响应放大部件的组件。这有效地提高了高频噪声的水平。

hfvt = fvtool(df，[1 -1]，1，'MagnitudeDisplay',“零”,'FS'Fs);传奇(hfvt“50阶FIR微分”,“DIFF函数的响应”);

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

V1 = DIFF（漂移）/ dt的;A1 = DIFF（V1）/ dt的;V1 = [0;V1];A1 = [0;0;A1];

区分使用第50阶FIR滤波器和补偿延迟。

D =意味着(grpdelay (df));%过滤器延迟V2 =滤波器（DF，[漂移;零（d，1）]）;V2 = V2（d + 1：结束）;A2 =滤波器（DF，[V2;零（d，1）]）;A2 = A2（d + 1：结束）;V2 = V2 / dt的;A2 = A2 / dt的^ 2;

绘制几个楼层位移数据点。此外，还绘制了一些数据点的速度和加速度计算与diff和50阶FIR滤波器。注意噪声是如何在速度估计值中被轻微放大，而在加速度估计值中被放大的diff。

helperFilterIntroductionPlot2（吨，漂移，V1，V2，A1，A2）

集成的信号

泄漏积分滤波器是具有传递函数的全极点滤波器 $H（Z）= 1 / [1-CZ ^ { - 1}]$ 在哪里 $ C $ 是一个常数，必须小于1才能保证滤波器的稳定性。这并不奇怪 $ C $ 接近1，则漏积分接近的逆diff转换功能。应用泄漏积分上一节分别取回速度和漂移中获得的加速度和速度估计。使用与所获得的估计diff功能，因为它们是噪音。

使用漏电积分器美元= 0.999 $ 。画出泄漏积分器滤波器的幅值响应。注意，该滤波器作为一个低通滤波器有效地消除高频噪声。

关闭（hfvt）fvtool（1，[1 -.999]'FS'，FS）

过滤速度和加速度与漏积分。

v_original = V1;a_original = A1;d_leakyint =过滤器（1，[1 -0.999]，v_original）;v_leakyint =过滤器（1，[1 -0.999]，a_original）;通过时间微分％乘法d_leakyint = d_leakyint * dt的;v_leakyint = v_leakyint * dt的;

画出位移和速度估计和比较原始的信号V1和A1。

helperFilterIntroductionPlot3 (t,漂移,d_leakyint、v_original v_leakyint)

您还可以使用集成的信号cumsum和cumtrapz功能。结果将类似于与漏积分获得的那些。

结论

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