今天學(xué)習(xí)的是 Signal Processing Toolbox的 一個(gè) Signal Smoothing 例程
Step1 .
加載數(shù)據(jù)哑蔫，使用的是波士頓的溫度數(shù)據(jù)碰声，每1小時(shí)一個(gè)數(shù)據(jù)姑蓝，共計(jì)31天，744個(gè)數(shù)據(jù)嘶是，然后畫圖
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)');
結(jié)果如下
Step2.
做滑動(dòng)平均债蜜，因?yàn)楦信d趣的是每一天對(duì)整個(gè)溫度的影響晴埂。所以用一個(gè)滑動(dòng)平均做濾波處理究反， 滑動(dòng)平均濾波器的長度為24，即24小時(shí)的長度儒洛。
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)');
處理后的效果如下
Step3.
去除移動(dòng)延遲精耐。對(duì)于每個(gè)對(duì)稱的長度為N的濾波器，其時(shí)間延遲為（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)');
效果如下
Step4.
抽取平均差異卦停。 查看每個(gè)時(shí)刻的溫度對(duì)整體溫度的影響， 所以把31天相同時(shí)間的溫度做平均處理恼蓬。
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函數(shù)后面有矩陣轉(zhuǎn)置
Step5.
用帶權(quán)重的滑動(dòng)平均濾波器濾波。用的是二項(xiàng)式濾波滚秩，(1/2,1/2)^n 专执，n很大時(shí)，接近于正太曲線郁油，當(dāng)n小時(shí)可以用來濾除高頻分量本股。濾波器系數(shù)的生成方法為，(1/2,1/2)和(1/2,1/2)做卷積桐腌，然后再繼續(xù)和(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)');
濾波效果如下
得到的濾波器系數(shù)如下
Step 6.
指數(shù)滑動(dòng)濾波，和高斯擴(kuò)展濾波類似案站。這種帶權(quán)重的滑動(dòng)濾波器的特點(diǎn)就是需要的窗口很小躬审，對(duì)資源占用小
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)');
從如下的放大波形可以看到，極值被去掉了
Step7.
為了更緊密的跟蹤信號(hào)蟆盐，可以采用 擬合式的濾波器來濾波承边， 這也是一種最小方差概念意義的濾波器。
采用sgolayfilt 這個(gè)函數(shù)來實(shí)現(xiàn)一個(gè) Savitzky-Golay smoothing filter石挂， 這個(gè)函數(shù)會(huì)自動(dòng)的計(jì)算參數(shù)和時(shí)間延遲等博助。其中第三個(gè)參數(shù)必須為奇數(shù)
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]);
Step8.
關(guān)于重采樣。
有時(shí)候?yàn)榱藨?yīng)用一個(gè)滑動(dòng)濾波器痹愚，可以先對(duì)信號(hào)重新采樣富岳。
先加載一個(gè)電壓信號(hào)
load openloop60hertz
fs = 1000;
t = (0:numel(openLoopVoltage)-1) / fs;
plot(t,openLoopVoltage);
ylabel('Voltage (V)');
xlabel('Time (s)');
title('Open-loop Voltage Measurement');
可以看到，信號(hào)采樣1KHz拯腮，但是有很多60Hz的交流干擾窖式。
先應(yīng)用一個(gè)一致的滑動(dòng)平均濾波器，因?yàn)?1000 / 60 = 16.667动壤， 所以滑動(dòng)窗大小取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');
濾波后效果如下
可以看到60Hz的干擾依然存在
因?yàn)?7 * 60 Hz = 1020 Hz萝喘，所以重采樣的頻率設(shè)為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');
可以看出干擾基本沒有了。
Step8.
中值濾波蜒灰。 有時(shí)候信號(hào)里面有毛刺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)');
重采樣并采用 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');
可見毛刺沒有去掉，且變寬了强窖。
采用中值濾波，因?yàn)槊痰膶挾戎挥?削祈，所以中值濾波的寬度設(shè)為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');
效果如下
總結(jié)：
學(xué)習(xí)到的函數(shù)主要有如下幾個(gè)，
filter髓抑，用來做濾波
reshape咙崎，重新排列數(shù)據(jù)
conv，卷積函數(shù)
sgolayfilt吨拍，用于 Savitzky-Golay smoothing filter
numel褪猛， 數(shù)組元素個(gè)數(shù)
medfilt1， 中值濾波幾種滑動(dòng)濾波方法
滑動(dòng)平均濾波羹饰，用來濾除高頻
二項(xiàng)式濾波伊滋，當(dāng)n比較小時(shí)，濾除高頻成分
指數(shù)式移動(dòng)平均濾波队秩，最濾波窗的要求小
Savitzky-Golay smoothing filter笑旺，多項(xiàng)式濾波
中值濾波，用來濾除毛刺為了能更好的應(yīng)用滑動(dòng)平均濾波馍资，可以先對(duì)數(shù)據(jù)進(jìn)行重采樣筒主。
滑動(dòng)平均濾波后的信號(hào)都有延遲，要把延遲考慮進(jìn)去鸟蟹。
