一维信号去噪分析.doc

一维信号去噪分析 一维信号滤波分析读取实验数据，绘制一维时域谱，如图1所示 图1下面对信号进行FFT变换：Fs=1024;%假定采样频率N=597;%实际采样点数y=xlsread(C:Documents and SettingsAdministrator桌面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/Nfigureplot(F(1:N/2),Ay2(1:N/2); %显示换算后的FFT模值结果title(幅度-频率曲线图);得到信号频谱如图2 图2 由于大多数信号噪声频率都为高频，所以可以先尝试设计一个简单的带通滤波器，这里采用窗函数法，根据FFT变换得到的频谱图取通带为50350Hz，算法函数如下：Fs=1024;N=597;n=36;y=xlsread(C:Documents and SettingsAdministrator桌面Book2.xls);y=y;x=xlsread(C:Documents and SettingsAdministrator桌面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:Documents and SettingsAdministrator桌面Book2.xls);y=y;x=xlsread(C:Documents and SettingsAdministrator桌面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=ca1 cdd1;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=0 0.12 0.13 0.20 0.21 0.58 0.59 0.65 0.66 1;m=0 0 1 1 0 0 1 1 0 0;b a=firls(n,f,m);滤波器频率响应如图7 图7先对用该滤波器对信号进行初步滤波，然后用平滑函数进行处理，最后用前面的小波分析得到滤波信号，其算法如下：y=xlsread(C:Documents and SettingsAdministrator桌面Book2.xls);y=y;x=xlsread(C:Documents and SettingsAdministrator桌面Book1.xls);x=xn=90;f=0 0.12 0.13 0.20 0.21 0.58 0.59 0.65 0.66 1;m=0 0 1 1 0 0 1 1 0 0;b a=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:Documents and SettingsAdministrator桌面Book2.xls);yf=y;x=xlsread(C:Documents and SettingsAdministrator桌面Book1.xls);x=x;yf_in=yf zeros(1,delay);%当前的输入信号yf_delay=zeros(1,delay),yf;%延迟的输入信号M=num+delay;y2=zeros(1,M);%输出初始值w=zeros(1,N);%滤波系数初始值for n=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;%系数调整endfigure;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 图9clear allclose allclcV1 = rand(20,1)*700+200;subplot 231; plot(V1); ylim(200 700); grid;y=mean5_3(V1,5);subplot 232; plot(y); ylim(200 700); grid;y=mean5_3(V1,10);subplot 233; plot(y); ylim(000 700); grid;y=mean5_3(V1,50);subplot 234; plot(y); ylim(200 700); grid;y=mean5_3(V1,200);subplot 235; plot(y); ylim(200 700); grid;y=mean5_3(V1,500);subplot 236; plot(y); ylim(200 700); grid;%-function y = mean5_3(x, m)% x为被处理的数据% m 为循环次数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;endy =a;subplot 231;plot(V1);grid;y=mean5_3(V1,5);subplot 232; plot(y);grid;y=mean5_3(V1,10);subplot 233;plot(y);grid;y=mean5_3(V1,50);subplot 234; plot(y);grid;y=mean5_3(V1,200);subplot 235;plot(y);grid;y=mean5_3(V1,500);subplot 236; plot(y);grid;