MATLAB学习help之——Signal Smoothing from Signal Processing Toolbox

今天学习的是 Signal Processing Toolbox的 一个 Signal Smoothing 例程

Step1 .
加载数据,使用的是波士顿的温度数据,每1小时一个数据,共计31天,744个数据,然后画图

load bostemp
days = (1:31*24)/24;
plot(days, tempC), axis tight;
ylabel('Temp (\circC)');
xlabel('Time elapsed from Jan 1, 2011 (days)');
title('Logan Airport Dry Bulb Temperature (source: NOAA)');

结果如下


图片.png

Step2.
做滑动平均,因为感兴趣的是每一天对整个温度的影响。所以用一个滑动平均做滤波处理, 滑动平均滤波器的长度为24,即24小时的长度。

hoursPerDay = 24;
coeff24hMA = ones(1, hoursPerDay)/hoursPerDay;

avg24hTempC = filter(coeff24hMA, 1, tempC);
plot(days, [tempC avg24hTempC]);
legend('Hourly Temp','24 Hour Average (delayed)','location','best');
ylabel('Temp (\circC)');
xlabel('Time elapsed from Jan 1, 2011 (days)');
title('Logan Airport Dry Bulb Temperature (source: NOAA)');

处理后的效果如下


图片.png

Step3.
去除移动延迟。对于每个对称的长度为N的滤波器,其时间延迟为(N-1)/2,所以

fDelay = (length(coeff24hMA)-1)/2;
plot(days, tempC, 'b', ...
     days-fDelay/24, avg24hTempC, 'g');
axis tight;
legend('Hourly Temp','24 Hour Average','location','best');
ylabel('Temp (\circC)');
xlabel('Time elapsed from Jan 1, 2011 (days)');
title('Logan Airport Dry Bulb Temperature (source: NOAA)');

效果如下


图片.png

Step4.
抽取平均差异。 查看每个时刻的温度对整体温度的影响, 所以把31天相同时间的温度做平均处理。

figure
deltaTempC = tempC - avg24hTempC;
deltaTempC = reshape(deltaTempC, 24, 31).';

plot(1:24, mean(deltaTempC)), axis tight;
title('Mean temperature differential from 24 hour average');
xlabel('Hour of day (since midnight)');
ylabel('Temperature difference (\circC)');

效果如下,注意reshape函数后面有矩阵转置


图片.png

Step5.
用带权重的滑动平均滤波器滤波。用的是二项式滤波,(1/2,1/2)^n ,n很大时,接近于正太曲线,当n小时可以用来滤除高频分量。滤波器系数的生成方法为,(1/2,1/2)和(1/2,1/2)做卷积,然后再继续和(1/2,1/2)递归卷积。如下

h = [1/2 1/2];
binomialCoeff = conv(h,h);
for n = 1:4
    binomialCoeff = conv(binomialCoeff,h);
end

figure
fDelay = (length(binomialCoeff)-1)/2;
binomialMA = filter(binomialCoeff, 1, tempC);
plot(days, tempC, 'b', ...
     days-fDelay/24,  binomialMA, 'g');
axis tight;
legend('Hourly Temp','Binomial Weighted Average','location','best');
ylabel('Temp (\circC)');
xlabel('Time elapsed from Jan 1, 2011 (days)');
title('Logan Airport Dry Bulb Temperature (source: NOAA)');

滤波效果如下


图片.png

得到的滤波器系数如下


图片.png

Step 6.
指数滑动滤波,和高斯扩展滤波类似。这种带权重的滑动滤波器的特点就是需要的窗口很小,对资源占用小

alpha = 0.45;
exponentialMA = filter(alpha, [1 alpha-1], tempC);
plot(days, tempC, 'b', ...
     days-fDelay/24,  binomialMA, 'g', ...
     days-1/24,       exponentialMA,'r');

axis tight;
legend('Hourly Temp', ...
       'Binomial Weighted Average', ...
       'Exponential Weighted Average','location','best');
ylabel('Temp (\circC)');
xlabel('Time elapsed from Jan 1, 2011 (days)');
title('Logan Airport Dry Bulb Temperature (source: NOAA)');
图片.png

从如下的放大波形可以看到,极值被去掉了


图片.png

Step7.
为了更紧密的跟踪信号,可以采用 拟合式的滤波器来滤波, 这也是一种最小方差概念意义的滤波器。
采用sgolayfilt 这个函数来实现一个 Savitzky-Golay smoothing filter, 这个函数会自动的计算参数和时间延迟等。其中第三个参数必须为奇数

cubicMA   = sgolayfilt(tempC, 3, 7);
quarticMA = sgolayfilt(tempC, 4, 7);
quinticMA = sgolayfilt(tempC, 5, 9);
plot(days, [tempC cubicMA quarticMA quinticMA]);
legend('Hourly Temp','Cubic-Weighted MA', 'Quartic-Weighted MA', ...
       'Quintic-Weighted MA','location','southeast');
ylabel('Temp (\circC)');
xlabel('Time elapsed from Jan 1, 2011 (days)');
title('Logan Airport Dry Bulb Temperature (source: NOAA)');
axis([3 5 -5 2]);
图片.png

Step8.
关于重采样。
有时候为了应用一个滑动滤波器,可以先对信号重新采样。
先加载一个电压信号

load openloop60hertz
fs = 1000;
t = (0:numel(openLoopVoltage)-1) / fs;
plot(t,openLoopVoltage);
ylabel('Voltage (V)');
xlabel('Time (s)');
title('Open-loop Voltage Measurement');
图片.png

可以看到,信号采样1KHz,但是有很多60Hz的交流干扰。
先应用一个一致的滑动平均滤波器,因为 1000 / 60 = 16.667, 所以滑动窗大小取17

plot(t,sgolayfilt(openLoopVoltage,1,17));
ylabel('Voltage (V)');
xlabel('Time (s)');
title('Open-loop Voltage Measurement');
legend('Moving average filter operating at 58.82 Hz', ...
       'location','southeast');

滤波后效果如下


图片.png

可以看到60Hz的干扰依然存在

因为17 * 60 Hz = 1020 Hz,所以重采样的频率设为1020,如下

fsResamp = 1020;
vResamp = resample(openLoopVoltage, fsResamp, fs);
tResamp = (0:numel(vResamp)-1) / fsResamp;
vAvgResamp = sgolayfilt(vResamp,1,17);
plot(tResamp,vAvgResamp);
ylabel('Voltage (V)');
xlabel('Time (s)');
title('Open-loop Voltage Measurement');
legend('Moving average filter operating at 60 Hz','location','southeast');
图片.png

可以看出干扰基本没有了。

Step8.
中值滤波。 有时候信号里面有毛刺spikes,如下

spikeSignal = zeros(size(openLoopVoltage));
spikeSignal(1:100:2000) = -6;
noisyLoopVoltage = openLoopVoltage + spikeSignal;
plot(t, noisyLoopVoltage)
ylabel('Voltage (V)');
xlabel('Time (s)');
title('Open-loop Voltage Measurement (spikes added)');
图片.png

重采样并采用 Savitzky-Golay 滤波后的效果如下

vResamp = resample(noisyLoopVoltage, fsResamp, fs);
tResamp = (0:numel(vResamp)-1) / fsResamp;
vAvgResamp = sgolayfilt(vResamp,1,17);
plot(tResamp,vAvgResamp);
ylabel('Voltage (V)');
xlabel('Time (s)');
title('Open-loop Noisy Voltage Measurement (spikes added)');
legend('Moving average filter operating at 60 Hz','location','southeast');
图片.png

可见毛刺没有去掉,且变宽了。

采用中值滤波,因为毛刺的宽度只有1,所以中值滤波的宽度设为3,如下

medfiltLoopVoltage = medfilt1(noisyLoopVoltage, 3);
vResamp = resample(medfiltLoopVoltage, fsResamp, fs);
tResamp = (0:numel(vResamp)-1) / fsResamp;
vAvgResamp = sgolayfilt(vResamp,1,17);
plot(tResamp,vAvgResamp);
ylabel('Voltage (V)');
xlabel('Time (s)');
title('Open-loop Noisy Voltage Measurement (spikes added)');
legend('Median and moving average filtered','location','southeast');

效果如下


图片.png

总结:

  1. 学习到的函数主要有如下几个,
    filter,用来做滤波
    reshape,重新排列数据
    conv,卷积函数
    sgolayfilt,用于 Savitzky-Golay smoothing filter
    numel, 数组元素个数
    medfilt1, 中值滤波

  2. 几种滑动滤波方法
    滑动平均滤波,用来滤除高频
    二项式滤波,当n比较小时,滤除高频成分
    指数式移动平均滤波,最滤波窗的要求小
    Savitzky-Golay smoothing filter,多项式滤波
    中值滤波,用来滤除毛刺

  3. 为了能更好的应用滑动平均滤波,可以先对数据进行重采样。

  4. 滑动平均滤波后的信号都有延迟,要把延迟考虑进去。

推荐阅读更多精彩内容