一维信号去噪解析.docx
《一维信号去噪解析.docx》由会员分享，可在线阅读，更多相关《一维信号去噪解析.docx（11页珍藏版）》请在冰豆网上搜索。
一维信号去噪解析
一维信号滤波分析
读取实验数据，绘制一维时域谱，如图1所示
图1
下面对信号进行FFT变换：
Fs=1024;%假定采样频率
N=597;%实际采样点数
y=xlsread（'C:
\DocumentsandSettings\Administrator\桌面\Book2.xls'）;
y=y';
y1=fft（y）;
y2=10*log10（（abs（y1）.^2））;
figure;
plot（y2）;
title（'fft变换'）%做FFT变换
Ay2=abs（y2）;%取模
Ay2=Ay2/（N/2）;%换算成实际的幅度
Ay2
（1）=Ay2
（1）/2;
F=（[1:
N]-1）*Fs/N;%换算成实际的频率值，Fn=（n-1）*Fs/N
figure
plot（F（1:
N/2）,Ay2（1:
N/2））;%显示换算后的FFT模值结果
title（'幅度-频率曲线图'）;
得到信号频谱如图2
图2
由于大多数信号噪声频率都为高频，所以可以先尝试设计一个简单的带通滤波器，这里采用窗函数法，根据FFT变换得到的频谱图取通带为50~350Hz，算法函数如下：
Fs=1024;
N=597;
n=36;
y=xlsread（'C:
\DocumentsandSettings\Administrator\桌面\Book2.xls'）;
y=y';
x=xlsread（'C:
\DocumentsandSettings\Administrator\桌面\Book1.xls'）;
x=x';
Wn=[50,350]*2/Fs;%归一化频率
b=fir1（n,Wn）;
yf=filter（b,1,y）;
figure;
plot（x,y）;
xlabel（'运动位置'）;
ylabel（'信号幅值'）;
title（'窗函数法滤波后的信号'）
滤波器频谱响应如图3
滤波后得到滤波信号如图4
图4
从图4可以看出，噪声部分被滤除，但这样设计的带通滤波器对去噪没有明显的效果，并且在未知采样频率的情况下对不确定信号进行FFT变换的效果并不理想，而小波分析则克服了这种缺陷，它对随机信号具有多分辨率分析的特点，在时域和频域都有表征信号局部信息的能力，时间窗和频率窗都可以根据信号的具体形态动态调整，在一般情况下，在低频部分可以采用较低的时间分辨率，而提高频率的分辨率，在高频情况下可以用较低的频率分辨率来换取精确的时间定位。
针对以上实验数据，这里采用的思路是先对信号进行低层分解强度去噪，然后通过阈值信号处理，得到最终的滤波信号，小波算法函数设计如下：
y=xlsread（'C:
\DocumentsandSettings\Administrator\桌面\Book2.xls'）;
y=y';
x=xlsread（'C:
\DocumentsandSettings\Administrator\桌面\Book1.xls'）;
x=x';
[c,l]=wavedec（y,1,'db1'）;%用db1小波对原始信号进行1层分解并提取系数
ca1=appcoef（c,l,'db1',1）;%提取近似系数
cd1=detcoef（c,l,1）;%提取1阶详细系数
cdd1=zeros（1,length（cd1））;%对信号进行强制性去噪处理并显示结果
c1=[ca1cdd1];
s1=waverec（c1,l,'db1'）;%多层次的一维小波重构
N=2;%N越大，分解的层数越多，则能滤除更多噪声，同时取值过大也会导致信号失真
yd=wden（s1,'rigrsure','h','mln',N,'sym20'）;%rigrsure阈值信号处理
figure;
plot（x,yd）;
xlabel（'运动位置'）;
ylabel（'信号幅值'）;
title（'小波分析滤波后的信号'）
得到滤波信号如图5
图5
如图5所示，信号有明显改善，大部分噪声已去除，但是波形还是很不规则。
假设信号是由多个单频信号组成的组合信号，我们可以针对信号的频谱直接设计一个滤波器，而通过对信号进行FFT变换并不能得到想要的结果，所以这里采用另外一种谱分析方法——功率谱估计，它是专门针对随机信号所采取的谱分析方法。
功率谱估计分为非参数方法、参数方法、子空间方法三大类，子空间方法相对来说有更高的分辨率，它又可分为MUSIC方法和特征向量法。
其中，特征向量法主要适用于混有高斯白噪声的正弦信号的功率谱估计，而MUSIC法适合更为普遍情况下正弦信号的功率谱估计，所以这里采用MUSIC法对信号进行初步分析，得到功率谱图如图6
图6
由图6可以看出，信号在频率为0.16，0.62两点处分别出现峰值，我们假设信号就是由这两种频率构成的单频信号的组合，针对此信号频率设计一个滤波器，其基本步骤如下：
n=50;
f=[00.120.130.200.210.580.590.650.661];
m=[0011001100];
[ba]=firls（n,f,m）;
滤波器频率响应如图7
图7
先对用该滤波器对信号进行初步滤波，然后用平滑函数进行处理，最后用前面的小波分析得到滤波信号，其算法如下：
y=xlsread（'C:
\DocumentsandSettings\Administrator\桌面\Book2.xls'）;
y=y';
x=xlsread（'C:
\DocumentsandSettings\Administrator\桌面\Book1.xls'）;
x=x'
n=90;
f=[00.120.130.200.210.580.590.650.661];
m=[0011001100];
[ba]=firls（n,f,m）;
y1=filter（b,a,y）;
y1（1,[1:
100]）=y1（1,[371:
470]）;
y2=smooth（y1）;
y3=wden（y2,'rigrsure','h','mln',2,'sym20'）;
figure;
plot（x,y3）;
xlabel（'运动位置'）;
ylabel（'信号幅值'）;
title（'滤波信号'）
得到的滤波信号如图8
图8
由此得到的的滤波信号虽然规则平滑，可是信号幅值与实验数据相比有较大差距，所以此种方法的可行性有待进一步研究
此外，也可以采用LMS算法进行自适应滤波，得到的滤波信号也比较理想，其算法步骤如下：
num=597;%采样点数
delay=10;%延迟间隔
N=10;%滤波器抽头数
u=0.001;%调整步长
k=1:
num;
y=xlsread（'C:
\DocumentsandSettings\Administrator\桌面\Book2.xls'）;
yf=y';
x=xlsread（'C:
\DocumentsandSettings\Administrator\桌面\Book1.xls'）;
x=x';
yf_in=[yfzeros（1,delay）];%当前的输入信号
yf_delay=[zeros（1,delay）,yf];%延迟的输入信号
M=num+delay;
y2=zeros（1,M）;%输出初始值
w=zeros（1,N）;%滤波系数初始值
forn=N:
M-N+1
yfyf=yf_delay（n+N-1:
-1:
n）;%滤波器输入值
y2（n）=w*yfyf';%滤波器输出值
e（n）=yf_in（n）;
w=w+2*u.*e（n）.*yfyf;%系数调整
end
figure;
y3=e（10:
num）;%提取信号
y4=wden（y3,'rigrsure','h','mln',2,'sym20'）;%小波进一步滤波
y4=smooth（y4）;
x=x（:
[1:
588]）;%保持向量维数一致
plot（x,y4）;
xlabel（'运动位置'）;
ylabel（'信号幅值'）;
title（'LMS算法滤波后的信号'）;
得到的处理结果如图9
图9
clearall
closeall
clc
V1=rand（20,1）*700+200;
subplot231;plot（V1）;ylim（[200700]）;grid;
y=mean5_3（V1,5）;
subplot232;plot（y）;ylim（[200700]）;grid;
y=mean5_3（V1,10）;
subplot233;plot（y）;ylim（[000700]）;grid;
y=mean5_3（V1,50）;
subplot234;plot（y）;ylim（[200700]）;grid;
y=mean5_3（V1,200）;
subplot235;plot（y）;ylim（[200700]）;grid;
y=mean5_3（V1,500）;
subplot236;plot（y）;ylim（[200700]）;grid;
%------------------------------------------------------------
functiony=mean5_3（x,m）
%x为被处理的数据
%m为循环次数
n=length（x）;
a=x;
fork=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;
forj=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;
subplot231;
plot（V1）;
grid;
y=mean5_3（V1,5）;
subplot232;plot（y）;
grid;
y=mean5_3（V1,10）;
subplot233;
plot（y）;
grid;
y=mean5_3（V1,50）;
subplot234;plot（y）;
grid;
y=mean5_3（V1,200）;
subplot235;
plot（y）;
grid;
y=mean5_3（V1,500）;
subplot236;plot（y）;
grid;
升级会员 