信號(hào)處理實(shí)驗(yàn)-FIR維納濾波器

實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)題目假設(shè)一個(gè)點(diǎn)目標(biāo)在x，y平面上繞單位圓做圓周運(yùn)動(dòng)，由于外界干擾，其運(yùn)動(dòng)軌跡發(fā)生了偏移。其中，x方向的干擾為均值為0，方差為0.05的高斯噪聲；y方向干擾為均值為0，方差為0.06的高斯噪聲。產(chǎn)生滿足要求的x方向和y方向隨機(jī)噪聲500個(gè)樣本；明確期望信號(hào)和觀測(cè)信號(hào)；試設(shè)計(jì)一FIR維納濾波器，確定最佳傳遞函數(shù)：，并用該濾波器處理觀測(cè)信號(hào)，得到其最佳估計(jì)。（注：自行設(shè)定誤差判定閾值，根據(jù)閾值確定濾波器的階數(shù)或傳遞函數(shù)的長(zhǎng)度）。要求3中，也可以選擇Kalman濾波器進(jìn)行濾波處理，采用哪種濾波器可以自由選擇。分別繪制出x方向和y方向的期望信號(hào)、噪聲信號(hào)、觀測(cè)信號(hào)、濾波后信號(hào)、誤差信號(hào)的曲線圖；在同一幅圖中繪制出期望信號(hào)、觀測(cè)信號(hào)和濾波后點(diǎn)目標(biāo)的運(yùn)動(dòng)軌跡。實(shí)驗(yàn)報(bào)告要求：給出求解思路和結(jié)果分析，給出MATLAB實(shí)現(xiàn)源程序和程序注解。解題思路目標(biāo)信號(hào)在x、y方向均受到高斯噪聲的干擾，本實(shí)驗(yàn)分別從兩個(gè)方向分別進(jìn)行濾波處理，最終得到整體濾波處理后的信號(hào)。具體解題思路如下：選用FIR維納濾波器進(jìn)行濾波處理產(chǎn)生x、y方向噪聲信號(hào)xnoise和ynoise，產(chǎn)生x、y方向期望信號(hào)xs、ys，產(chǎn)生x、y方向觀測(cè)信號(hào)x、y。分別求出x、y方向上觀測(cè)信號(hào)的自相關(guān)矩陣以及觀測(cè)信號(hào)與期望信號(hào)的互相關(guān)矩陣，根據(jù)公式確定最佳傳遞函數(shù)。根據(jù)=分別求出x、y方向上最小均方誤差。MATLAB實(shí)現(xiàn)源程序closeall;clearall;%x軸方向和y軸方向期望信號(hào)x0=0;y0=0;r=1;theta=0:pi/249.5:2*pi;xs=x0+r*cos(theta);ys=y0+r*sin(theta);%產(chǎn)生x軸方向和y軸方向噪聲n=500;%產(chǎn)生500個(gè)隨機(jī)樣本ux=0;vx=0.05;t=randn(1,n);%產(chǎn)生1*n個(gè)偽隨機(jī)數(shù)t=t/std(t);t=t-mean(t);xnoise=ux+sqrt(vx)*t;uy=0;vy=0.06;g=randn(1,n);g=g/std(g);g=g-mean(g);ynoise=uy+sqrt(vy)*g;%x軸方向和y軸方向觀測(cè)信號(hào)x=xs+xnoise;y=ys+ynoise;%x方向上觀測(cè)信號(hào)的自相關(guān)矩陣rxx=xcorr(x);fori=1:100forj=1:100mrxx(i,j)=rxx(500-i+j);endend%x方向上觀測(cè)信號(hào)與期望信號(hào)的互相關(guān)矩陣rxs=xcorr(x,xs);fori=1:100mrxs(i)=rxs(499+i);endhoptx=inv(mrxx)*mrxs';%求出x方向上最佳傳遞函數(shù)fx=conv(x,hoptx);%維納濾波后x方向上的輸出nx=sum(abs(xs).^2);eminx=nx-mrxs*hoptx;%求出x方向上的最小均方誤差%y方向上觀測(cè)信號(hào)的自相關(guān)矩陣ryy=xcorr(y);fori=1:100forj=1:100mryy(i,j)=ryy(500-i+j);endend%y方向上觀測(cè)信號(hào)與期望信號(hào)的互相關(guān)矩陣rys=xcorr(y,ys);fori=1:100mrys(i)=rys(499+i);endhopty=inv(mryy)*mrys';%求出y方向上最佳傳遞函數(shù)fy=conv(y,hopty);%維納濾波后y方向上的輸出ny=sum(abs(ys).^2);eminy=ny-mrys*hopty;%求出y方向上的最小均方誤差plot(xs,ys,'-r',x0,y0,'.');holdon;plot(x,y,'b:');holdon;plot(fx,fy,'g-');legend('期望信號(hào)','圓周運(yùn)動(dòng)中心','觀測(cè)信號(hào)','濾波后點(diǎn)目標(biāo)的運(yùn)動(dòng)軌跡')axissquare;figure;subplot(2,4,1)plot(xs);title('x方向期望信號(hào)');subplot(2,4,2)plot(xnoise);title('x方向噪聲信號(hào)');subplot(2,4,3)plot(x);title('x方向觀測(cè)信號(hào)');subplot(2,4,4)n=0:500;plot(n,eminx);title('x方向最小均方誤差');subplot(2,4,5)plot(ys);title('y方向期望信號(hào)');subplot(2,4,6)plot(ynoise);title('y方向噪聲信號(hào)');subplot(2,4,7)plot(y);title('y方向觀測(cè)信號(hào));subplot(2,4,8)plot(n,eminy);title('y方向最小均方誤差');結(jié)果及其分析運(yùn)行結(jié)果如下圖：分析：運(yùn)行結(jié)果中濾波后的到的信號(hào)與原始信號(hào)和噪聲信號(hào)的對(duì)比可以看出，濾波后的結(jié)果與期望信號(hào)還是很接近的，整體上達(dá)到了最優(yōu)濾波的效果。濾波后點(diǎn)目標(biāo)的運(yùn)動(dòng)軌跡有回到零點(diǎn)的情況發(fā)生，分析其原因是：傳遞函數(shù)的點(diǎn)數(shù)為M，x的點(diǎn)數(shù)為N，進(jìn)行卷積時(shí)，會(huì)自動(dòng)補(bǔ)0使二者長(zhǎng)度相同，卷積之后點(diǎn)數(shù)為Ｍ＋Ｎ－１。實(shí)驗(yàn)總結(jié)試驗(yàn)中遇到幾個(gè)問題，經(jīng)過查閱資料以及和同學(xué)討論最終得到解決。此過程中遇到的問題及解決方案如下：?jiǎn)栴}1：點(diǎn)目標(biāo)做圓周運(yùn)動(dòng)，將圓周中心定于非零點(diǎn)，經(jīng)過濾波后，均方誤差較大，濾波效果較差。解決：忽略了信號(hào)的零均值化處理，將圓周中心定于零點(diǎn)，這樣信號(hào)的均值為零，相當(dāng)于進(jìn)行了零均值化處理。問題2：由于產(chǎn)生x、y方向上隨機(jī)噪聲的方法相同，將上述程序中產(chǎn)生噪聲信號(hào)的部分改為如下：%產(chǎn)生x軸方向和y軸方向噪聲n=500;%產(chǎn)生500個(gè)隨機(jī)樣本ux=0;vx=0.05;uy=0;vy=0.06;t=randn(1,n);%產(chǎn)生1*n個(gè)偽隨機(jī)數(shù)t=t/std(t);t=t-mean(t);xnoise=ux+sqrt(vx)*t;ynoise=uy+sqrt(vy)*g;運(yùn)行結(jié)果如下圖，發(fā)現(xiàn)，x、y方向干擾加進(jìn)去后，造成的干擾有一定的規(guī)律，方向大概一致。解決：分析發(fā)現(xiàn)，產(chǎn)生的x、y方向上隨機(jī)噪聲用的是同一組偽隨機(jī)數(shù)，才出現(xiàn)下圖的結(jié)果。所以改為分別產(chǎn)生。問題3：圓周的等間隔采樣問題，采用的間隔為2*pi/50