Main Content

取信号的导数

您要在不增加噪声功率的情况下对信号求导。MATLAB® 提供的函数 diff 会放大噪声,对于高阶导数会恶化不精确性。要解决此问题,请改用微分滤波器。

分析地震时建筑物楼层的位移。找到速度和加速度作为时间的函数。

加载文件 earthquake。该文件包含以下变量:

  • drift:楼层位移,以厘米为单位进行测量

  • t:时间,以秒为单位进行测量

  • Fs:采样率,等于 1 kHz

load('earthquake.mat')

使用 pwelch 显示信号功率谱的估计值。请注意大部分信号能量包含在低于 100 Hz 的频率中。

pwelch(drift,[],[],[],Fs)

使用 designfilt 设计一个阶数为 50 的 FIR 微分器。要包含大部分信号能量,请指定 100 Hz 的通带频率和 120 Hz 的阻带频率。使用 zerophase 检查滤波器。

Nf = 50; 
Fpass = 100; 
Fstop = 120;

d = designfilt('differentiatorfir','FilterOrder',Nf, ...
    'PassbandFrequency',Fpass,'StopbandFrequency',Fstop, ...
    'SampleRate',Fs);

zerophase(d,[],Fs)

对漂移求导以求出速度。将导数除以 dt(即连续采样之间的时间间隔),以设置正确的单位。

dt = t(2)-t(1);

vdrift = filter(d,drift)/dt;

滤波后的信号存在延迟。使用 grpdelay 确定延迟是滤波器阶数的一半。通过丢弃采样对此进行补偿。

delay = mean(grpdelay(d))
delay = 25
tt = t(1:end-delay);
vd = vdrift;
vd(1:delay) = [];

输出还包括瞬变,其长度等于滤波器阶数,或者是群延迟的两倍。在上面已丢弃 delay 个采样。再次丢弃 delay 个采样以消除瞬变。

tt(1:delay) = [];
vd(1:delay) = [];

对漂移和漂移速度绘图。使用 findpeaks 验证漂移的最大值和最小值对应于其导数的过零点。

[pkp,lcp] = findpeaks(drift);
zcp = zeros(size(lcp));

[pkm,lcm] = findpeaks(-drift);
zcm = zeros(size(lcm));

subplot(2,1,1)
plot(t,drift,t([lcp lcm]),[pkp -pkm],'or')
xlabel('Time (s)')
ylabel('Displacement (cm)')
grid

subplot(2,1,2)
plot(tt,vd,t([lcp lcm]),[zcp zcm],'or')
xlabel('Time (s)')
ylabel('Speed (cm/s)')
grid

对漂移速度求微分以求出加速度。时滞长度是原来的两倍。丢弃两倍数量的采样来补偿延迟,丢弃相同数量的采样来消除瞬变。对速度和加速度绘图。

adrift = filter(d,vdrift)/dt;

at = t(1:end-2*delay);
ad = adrift;
ad(1:2*delay) = [];

at(1:2*delay) = [];
ad(1:2*delay) = [];

subplot(2,1,1)
plot(tt,vd)
xlabel('Time (s)')
ylabel('Speed (cm/s)')
grid

subplot(2,1,2)
plot(at,ad)
ax = gca;
ax.YLim = 2000*[-1 1];
xlabel('Time (s)')
ylabel('Acceleration (cm/s^2)')
grid

使用 diff 计算加速度。添加零来补偿数组大小的变化。将结果与使用滤波器获得的结果进行比较。请注意高频噪声的数量。

vdiff = diff([drift;0])/dt;
adiff = diff([vdiff;0])/dt;

subplot(2,1,1)
plot(at,ad)
ax = gca;
ax.YLim = 2000*[-1 1];
xlabel('Time (s)')
ylabel('Acceleration (cm/s^2)')
grid
legend('Filter')
title('Acceleration with Differentiation Filter')

subplot(2,1,2)
plot(t,adiff)
ax = gca;
ax.YLim = 2000*[-1 1];
xlabel('Time (s)')
ylabel('Acceleration (cm/s^2)')
grid
legend('diff')

另请参阅

| | | |

相关主题