從現(xiàn)實環(huán)境采集到的數(shù)據(jù)中經(jīng)常混疊有微弱噪聲，其中包括由于系統(tǒng)不穩(wěn)定產(chǎn)生的噪聲，也有周圍環(huán)境引入的毛刺，這些弱噪聲都需要在處理信號之前盡可能地消除或減弱。這一工作往往作為預處理的一部分。下面將介紹幾種簡單又實用的平滑處理方法：五點三次平滑法、MATLAB自帶平滑處理的smooth 函數(shù)和Savitzky-Golay平滑濾波器等。
一、五點三次平滑法
對于帶毛刺或弱噪聲的數(shù)據(jù)經(jīng)常會采用五點三次平滑法來進行平滑處理。
五點三次平滑法是利用最小二乘法原理對離散數(shù)據(jù)進行三次最小二乘多項式平滑的處理方法。
五點三次平滑法的函數(shù)為mean5_3:
函數(shù)：mean5_3
功能：對數(shù)據(jù)進行五點三次平滑處理
調(diào)用格式：
y=mean5_3(x,m)
說明：x是要平滑的輸入序列，m是對數(shù)據(jù)進行多次循環(huán)平滑的次數(shù)；y是平滑后的輸出序列。數(shù)據(jù)x能進行多次五點三次的平滑處理，但m必須選擇一個適當?shù)闹?，不宜太大，否則容易使峰值降低，峰值頻帶變寬。
函數(shù)程序如下：
function y=mean5_3(x,m)
% x為被處理的數(shù)據(jù)
% m 為循環(huán)次數(shù)
n=length(x);
a=x;
for k=1: m
b(1) = (69*a(1) +4*(a(2) +a(4)) -6*a(3) -a(5)) /70;
b(2) = (2* (a(1) +a(5)) +27*a(2) +12*a(3) -8*a(4)) /35;
for j=3:n-2
b (j) = (-3*(a(j-2) +a(j+2)) +12*(a(j-1) +a(j+1)) +17*a(j)) /35;
end
b (n-1) = (2*(a(n) +a(n-4)) +27*a(n-1) +12*a(n-2) -8*a(n-3)) /35;
b (n) = (69*a(n) +4* (a(n-1) +a(n-3)) -6*a(n-2) -a(n-4)) /70;
a=b;
end
y =a;案例1、有一組帶噪信號數(shù)據(jù)文件xnoisedata.txt,數(shù)據(jù)的第1列是時間,第2列是實驗檢測到的數(shù)據(jù)。由于環(huán)境的原因，實驗數(shù)據(jù)中含有噪聲，要求通過平滑方法對數(shù)據(jù)進行處理。程序如下：
clear all; clc; close all;
xx=load('xnoisedata1.txt'); % 讀入數(shù)據(jù)
time=xx(:,1); % 時間序列
x=xx(:,2); % 帶噪數(shù)據(jù)
xmean=mean5_3(x,50); % 調(diào)用mean5_3函數(shù),平滑數(shù)據(jù)
% 作圖
subplot 211; plot(time,x,'k');
xlabel('時間/s'); ylabel('幅值')
title('原始數(shù)據(jù)'); xlim([0 max(time)]);
subplot 212; plot(time,xmean,'k');
xlabel('時間/s'); ylabel('幅值')
title('平滑處理后的數(shù)據(jù)'); xlim([0 max(time)]);
set(gcf,'color','w');
運行結(jié)果如下：
?
二、MATLAB自帶的平滑函數(shù)smooth
MATLAB自帶的平滑函數(shù)smooth是一個低通濾波器，可見使用這種簡單形式的平均器，可以起到去除噪聲、提高信噪比的作用。這類濾波器被稱為平均濾波器或滑動平均濾波器。
MATALB自帶的平滑函數(shù)smooth
對MATLAB中自帶的平滑函數(shù)smooth介紹如下。
名稱：smooth
功能：對一維信號進行平滑處理
調(diào)用格式：
y=smooth(x)
y=smooth(x,SPAN)
y= smooth(x,SPAN,method)
說明：x是被平滑處理的輸入數(shù)據(jù)，SPAN是平滑處理中取的窗長（取奇數(shù)）。y是平滑處理后的輸出數(shù)據(jù)。Method是平滑處理的方法，共有6種，見下表。
案例2、有一個整周期的余弦信號，但疊加了噪聲，要求尋找余弦信號極小值的位置和幅值。調(diào)用MATLAB自帶的smooth函數(shù)試驗不同的參數(shù)，程序如下：
%
% pr4_5_2 from sy12
clear all; clc; close all;
k=1:500; % 產(chǎn)生一個從0到2*pi的向量，長為500
dn=2*pi/500;
x=cos((k-1)*dn); % 產(chǎn)生一個周期余弦信號
[val,loc]=min(x); % 求出余弦信號中的最小值幅值和位置
N=length(x); % 數(shù)據(jù)長
ns=randn(1,N); % 產(chǎn)生隨機信號
y=x+ns(1:N)*0.1; % 構(gòu)成帶噪信號
Err=var(x-y); % 求x-y的方差
fprintf('%4d %5.4f %5.6f\n',loc,val,Err);
y=y'; % 轉(zhuǎn)成列向量
z1=smooth(y,51,'moving'); % 'moving'平滑
Err1=var(x'-z1); % 求x-z1的方差
[v1,k1]=min(z1); % 求出平滑信號z1中的最小值幅值和位置
fprintf('1 %4d %5.4f %4d %5.6f\n',k1,v1,abs(loc-k1),Err1); % 顯示
z2=smooth(y,51,'lowess'); % 'lowess'平滑
Err2=var(x'-z2); % 求x-z2的方差
[v2,k2]=min(z2); % 求出平滑信號z2中的最小值幅值和位置
fprintf('2 %4d %5.4f %4d %5.6f\n',k2,v2,abs(loc-k2),Err2);
z3=smooth(y,51,'loess'); % 'loess'平滑
Err3=var(x'-z3); % 求x-z3的方差
[v3,k3]=min(z3); % 求出平滑信號z3中的最小值幅值和位置
fprintf('3 %4d %5.4f %4d %5.6f\n',k3,v3,abs(loc-k3),Err3);
z4=smooth(y,51,'sgolay',3); % 'sgolay'平滑
Err4=var(x'-z4); % 求x-z4的方差
[v4,k4]=min(z4); % 求出平滑信號z4中的最小值幅值和位置
fprintf('4 %4d %5.4f %4d %5.6f\n',k4,v4,abs(loc-k4),Err4);
z5=smooth(y,51,'rlowess'); % 'rlowess'平滑
Err5=var(x'-z5); % 求x-z5的方差
[v5,k5]=min(z5); % 求出平滑信號z5中的最小值幅值和位置
fprintf('5 %4d %5.4f %4d %5.6f\n',k5,v5,abs(loc-k5),Err5);
z6=smooth(y,51,'rloess'); % 'rloess'平滑
Err6=var(x'-z6); % 求x-z6的方差
[v6,k6]=min(z6); % 求出平滑信號z6中的最小值幅值和位置
fprintf('6 %4d %5.4f %4d %5.6f\n',k6,v6,abs(loc-k6),Err6);
% 作圖
subplot 211; plot(k,x,'k');
grid; xlim([0 500]); title('一周期余弦信號')
xlabel('樣點'); ylabel('幅值')
subplot 212; plot(k,y,'k'); %hold on
grid; axis([0 500 -1.5 1.5]); title('帶噪一周期余弦信號')
xlabel('樣點'); ylabel('幅值')
set(gcf,'color','w');
figure
subplot 321; plot(k,z1,'k'); title('moving法')
grid; axis([0 500 -1.5 1.5]); ylabel('幅值')
subplot 322; plot(k,z2,'k'); title('lowess法')
grid; axis([0 500 -1.5 1.5]); ylabel('幅值')
subplot 323; plot(k,z3,'k'); title('loess法')
grid; axis([0 500 -1.5 1.5]); ylabel('幅值')
subplot 324; plot(k,z4,'k'); title('sgolay法')
grid; axis([0 500 -1.5 1.5]); ylabel('幅值')
subplot 325; plot(k,z5,'k'); title('rlowess法')
grid; axis([0 500 -1.5 1.5]); xlabel('樣點'); ylabel('幅值')
subplot 326; plot(k,z6,'k'); title('rloess法')
grid; axis([0 500 -1.5 1.5]); xlabel('樣點'); ylabel('幅值')
set(gcf,'color','w');
運行結(jié)果如下：
?
?三、使用Savitzky-Golay函數(shù)平滑信號
使用Savitzky-Golay平滑濾波器對信號濾波，實際上是擬合了信號中的低頻成分，而將高頻成分“平滑”出去了。如果噪聲在高頻端，那么擬合的結(jié)果是去除了噪聲；反之，若噪聲在低頻端，信號在高頻端，那么濾波的結(jié)果是留下了噪聲。當然，用原信號減去噪聲，又可得到所期望的信號。
MATLAB中自帶的Savitzky-Golay平滑濾波函數(shù)
在MATLAB中自帶了兩個與Savitzky-Golay平滑濾波有關(guān)的函數(shù)，現(xiàn)介紹如下。
(1)、求 Savitzky-Golay 濾波器系數(shù)
名稱：sgolay
功能：設(shè)計低通濾波器求其系數(shù)
調(diào)用格式：
b= sgolay(k,f)
說明：輸入?yún)?shù)k是多項式擬合中的階數(shù)，f是窗長，該值必須為奇數(shù)。輸出參數(shù)b是Savitzky-Golay法的FIR平滑濾波器。
(2)、實現(xiàn) Savitzky-Golay 濾波
名稱：sgolayfilt
功能：實現(xiàn)Savitzky-Golay 濾波
調(diào)用格式：
y= sgolayfilt(x,k,f)
說明：輸入?yún)?shù)x是輸入信號，k是多項式擬合中的階數(shù)，f是窗長，該值必須為奇數(shù)。輸出參數(shù)y是FIR濾波器輸出。
案例3、 設(shè)正弦信號的采樣頻率為5Hz,信號頻率為0.2Hz,長200s,在信號中混入了噪聲，用Savitzky-Golay濾波器平滑該正弦信號，程序如下：
clear all; clc; close all;
t=0:.2:199; % 設(shè)置時間序列
s=10*sin(0.4*pi*t); % 原始信號
ns=randn(size(s)); % 產(chǎn)生噪聲序列
y=s+ns; % 構(gòu)成帶噪信號
x=sgolayfilt(y,3,19); % 通過Savitzky-Golay濾波器
% 作圖
figure
plot(t,y,'r');
xlim([0 20]); hold on; grid;
plot(t,x,'k');
xlabel('時間'); ylabel('幅值');
title('Savitzky-Golay濾波器的輸/輸出波形圖')
set(gcf,'color','w');運行結(jié)果如下：
利用Savitzky-Golay濾波器來消除噪聲時，一般窗長f值取得都不大，但Savitzky-Golay濾波器還可以求出信號的趨勢項，此時窗長f值取得都較大。
實驗數(shù)據(jù)下載鏈接如下：
https://mp.csdn.net/mp_download/manage/download/UpDetailed文章來源：http://www.zghlxwxcb.cn/news/detail-469494.html
參考文獻：MATLAB數(shù)字信號處理85個實用案例精講——入門到進階；宋知用（編著）文章來源地址http://www.zghlxwxcb.cn/news/detail-469494.html
到了這里，關(guān)于信號處理中簡單實用的方法——對信號進行平滑處理的文章就介紹完了。如果您還想了解更多內(nèi)容，請在右上角搜索TOY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！
