数字滤波实践介绍 -凯发k8网页登录

此示例说明如何设计、分析数字滤波器并将其应用于数据。它将帮助您回答以下问题：如何补偿滤波器引入的延迟？如何避免信号失真？如何从信号中去除不需要的内容？如对信号求微分？如何对信号求积分？

滤波器可用于以所需方式形成信号频谱，或执行数学运算，如微分和积分。接下来，您将了解在实践中使用的一些概念，以便在需要时轻松使用滤波器。

此示例重点介绍数字滤波器的应用，而不是其设计。如需了解有关如何设计数字滤波器的详细信息，请参阅示例。

补偿滤波引入的延迟

数字滤波器会在信号中引入延迟。根据滤波器特征，延迟可以在所有频率上保持不变，也可以随频率而变化。延迟的类型确定您必须采取的补偿措施。grpdelay 函数允许您以频率的函数形式查看滤波器延迟。通过查看此函数的输出，可以确定滤波器的延迟是常量还是随频率而变化（换句话说，是否与频率有关）。

通过按时间偏移信号，可以轻松补偿在所有频率上为常量的滤波器延迟。fir 滤波器通常具有常量延迟。另一方面，随频率变化的延迟会导致相位失真，并显著改变信号波形。与频率有关的延迟的补偿并不像常量延迟那样无关紧要。iir 滤波器会引入与频率有关的延迟。

补偿常量滤波器延迟

如前所述，您可以测量滤波器的群延迟，以验证它是频率的常量函数。您可以使用 grpdelay 函数来测量滤波器延迟 d，并通过在输入信号中追加 d 个零并按时间将输出信号偏移 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 个零。这可以确保所有有用的采样都通过滤波器，并且输入信号和经过延迟补偿的输出信号具有相同的长度。对数据进行滤波，并通过将输出信号偏移 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 个采样的等效滤波器对数据进行滤波。其他效应是，您得到的滤波器传递函数等于原始滤波器传递函数的平方幅值，并且滤波器阶数是原始滤波器阶数的两倍。

以上一节中定义的 ecg 信号为例。对此信号进行带延迟补偿和不带延迟补偿的滤波。设计一个截止频率为 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.

如果您的应用允许非因果前向/后向滤波运算，并且允许滤波器响应更改为原始响应的平方，则零相位滤波是很好的工具。

引入常时滞的滤波器是线性相位滤波器。引入与频率有关的延迟的滤波器是非线性相位滤波器。

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

滤波器通常用于从信号中去除不需要的频谱内容。您可以从各种滤波器中进行选择来实现以上目的。当您要去除高频成分时，可以选择低通滤波器；或者当您要去除低频成分时，可以选择高通滤波器。您也可以选择带通滤波器来去除低频和高频成分，同时保持中间频带不变。当您要去除给定频带上的频率时，可以选择带阻滤波器。

假设存在一个包含电力线嗡嗡声和白噪声的音频信号。电力线嗡嗡声是由 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 \times 1400 = 2800$ hz 就足以正确表示信号，但您使用的采样率为 44100 hz，这会浪费资源，因为您需要处理不必要的采样。您可以对信号进行下采样以降低采样率，并通过减少需要处理的采样数来减轻计算负载。较低的采样率还允许您设计更尖锐、更窄的带阻滤波器，以较小的滤波器阶数去除 60 hz 噪声。

按因子 10 对低通滤波信号进行下采样，以获得 fs/10 = 4.41 khz 采样率。绘制下采样之前与之后的信号频谱。

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 音调。将阻带宽度设置为 4 hz，以 60 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'})

如果您的计算机有声卡，可以聆听处理前后的信号。如上所述，最终结果是您有效地衰减了音频文件中的 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 函数能够对信号求微分，但其缺点是可能增加输出的噪声级别。更好的选择是使用微分滤波器，它在感兴趣的频带中充当微分器，在所有其他频率上充当衰减器，从而有效地去除高频噪声。

下面通过示例来分析地震过程中建筑物楼层的位移速度。在地震条件下，位移或漂移测量值记录在三层待测建筑物的第一楼层，并保存在 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

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

df = designfilt('differentiatorfir','filterorder',50,...
                'passbandfrequency',100,'stopbandfrequency',120,...
                'samplerate',fs);

diff 函数可被视为一阶 fir 滤波器，其响应为 $h (z) = 1 - z^{- 1}$ 。使用 fvtool 比较 50 阶 fir 微分滤波器的幅值响应和 diff 函数的响应。显然，这两种响应在通带区域（从 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 函数运算而缺失的采样。

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 ( c m / s squared baseline ) contains 2 objects of type line. these objects represent estimated acceleration using diff, estimated acceleration using fir filter.

对信号进行积分

泄漏积分器滤波器是一种全极点滤波器，它使用传递函数 $h (z) = 1 / [1 - c z^{- 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.

附录

此示例中使用以下辅助函数：

另请参阅

| | | | | | | |