今天学习的是 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
总结：
-
学习到的函数主要有如下几个，
filter，用来做滤波
reshape，重新排列数据
conv，卷积函数
sgolayfilt，用于 Savitzky-Golay smoothing filter
numel， 数组元素个数
medfilt1， 中值滤波 -
几种滑动滤波方法
滑动平均滤波，用来滤除高频
二项式滤波，当n比较小时，滤除高频成分
指数式移动平均滤波，最滤波窗的要求小
Savitzky-Golay smoothing filter，多项式滤波
中值滤波，用来滤除毛刺 -
为了能更好的应用滑动平均滤波，可以先对数据进行重采样。
-
滑动平均滤波后的信号都有延迟，要把延迟考虑进去。
网友评论