現(xiàn)代信號處置

4.信號的函數(shù)體現(xiàn)式為：，其中，為一隨時間變化的隨機過程，為通過390—410Hz帶通濾波器后的高斯白噪聲，為高斯白噪聲，采樣頻率為1kHz，采樣時間為2.048s。（1）運用當代信號解決的知識進行信號譜預計；（2）運用當代信號解決知識進行信號的頻率提?。唬?）分別運用Winner濾波和Kalman濾波進行去噪；（4）運用Wigner-Ville分布分析信號的時頻特性。（1）：運用當代信號解決的知識進行信號譜預計：典型譜預計中兩種重要的辦法為直接法和間接法，其中間接法則先根據(jù)N個樣本數(shù)據(jù)的樣本自有關函數(shù) （4.1）其中，且。計算樣本自有關函數(shù)的Fourier變換，得到功率譜 （4.2）周期圖辦法預計的功率譜為有偏預計，可通過加窗來減少其偏差。定義為 （4.3）式中 （4.4）式中，是窗函數(shù)的Fourier變換。功率譜預計程序為：clearclccloseallhiddensf=1000;nfft=2048;t=0:1/1000:2.047;A=normrnd(0,1,1,2048);N=wgn(1,2048,1);f1=390;f2=410;wc1=2*f1/sf;wc2=2*f2/sf;%歸一化頻率f0=[0wc1-0.05wc1wc2wc2+0.051];B=[001100];%設立帶通和帶阻weigh=[111];%設立帶通和帶阻權重b=remez(50,f0,B,weigh);%傳函分子D=filter(b,1,N);y=sin(2*pi*t*100)+1.5*sin(2*pi*t*300)+A.*sin(2*pi*t*200)+D+N;a(1,:)=y;a(2,:)=y.*sin(y);x=a(1,:);y=a(2,:)-a(1,:);f=0:sf/nfft:sf/2-sf/nfft;w=boxcar(nfft);%加矩形窗z=psd(y,nfft,sf,w,nfft/2);nn=1:nfft/2;plot(f(nn),abs(z(nn)));xlabel('頻率(Hz)');ylabel('幅值');gridon;圖4.1功率譜預計成果圖(2).信號頻率的提取用離散傅立葉算法離散傅立葉算法程序clearclccloseallhiddensf=1000;nfft=2048;t=0:1/1000:2.047;A=normrnd(0,1,1,2048);N=wgn(1,2048,1);f1=390;f2=410;wc1=2*f1/sf;wc2=2*f2/sf;%歸一化頻率f0=[0wc1-0.05wc1wc2wc2+0.051];B=[001100];%設立帶通和帶阻weigh=[111];%設立帶通和帶阻權重b=remez(50,f0,B,weigh);%傳函分子D=filter(b,1,N);y=sin(2*pi*t*100)+1.5*sin(2*pi*t*300)+A.*sin(2*pi*t*200)+D+N;t2=(0:nfft-1)/sf;f=(0:nfft-1)*sf/nfft;y1=abs(fft(y));f=f(1:nfft/2);y1=y1(1:nfft/2);plot(t,y);title('原始信號');axis([02.047-68]);plot(f,y1);title('fft頻率提取');axis([050001000]);xlabel('f/Hz');gridon;圖4.2原始信號時域圖圖4.3信號頻譜(3)分別運用Winner濾波和Kalman濾波進行去噪；clearallcloseallM=100;%維納濾波器階數(shù)sf=1000;nfft=2048;L=nfft;t=0:1/1000:2.047;A=normrnd(0,1,1,2048);N=wgn(1,2048,1);f1=390;f2=410;wc1=2*f1/sf;wc2=2*f2/sf;%歸一化頻率f0f0=[0wc1-0.05wc1wc2wc2+0.051];B=[001100];%設立帶通和帶阻weigh=[111];%設立帶通和帶阻權重b=remez(50,f0,B,weigh);%傳函分子D=filter(b,1,N);y=sin(2*pi*t*100)+1.5*sin(2*pi*t*300)+A.*sin(2*pi*t*200)+D+N;phixx=xcorr(y,y);fori=1:Mforj=1:MRxx(i,j)=phixx(i-j+L);endends=sin(2*pi*t*100)+1.5*sin(2*pi*t*300)+A.*sin(2*pi*t*200);phixs=xcorr(y,s);fori=1:Mrxs(i)=phixs(i+L);endh1=(inv(Rxx))*rxs';%獲得抱負FIR濾波器系數(shù)h1AA=sin(2*pi*t*100)+1.5*sin(2*pi*t*300)+A.*sin(2*pi*t*200);fori=1:Mh(i)=AA(i);end%繪圖比較預計濾波器與實際濾波器figurek=1:M;plot(k,h(k),'r',k,h1(k),'b');title('Idealh(n)&Calculatedh(n)');legend('Idealh(n)','Calculatedh(n)');xlabel('n');ylabel('h(n)');%比較抱負輸出與實際輸出v=D+N;S=conv(h,v);SI(1)=S(1);LL1=sin(2*pi*t*100)+1.5*sin(2*pi*t*300)+A.*sin(2*pi*t*200);fori=2:LSI(i)=LL1(i);endfigurek=1:L;plot(k,s(k),'r',k,SI(k),'b');title('s(n)VS.SI(n)');legend('s(n)','SI(n)',0);xlabel('n');ylabel('IdealOutput');holdonSR=conv(h1,y);figurek=1:L;plot(k,s(k),'r',k,SR(k),'b');title('s(n)VS.SR(n)');legend('s(n)去噪前','SR(n)去噪后',0);xlabel('n');ylabel('ActualOutput');圖4.4Winner濾波去噪圖Kalman濾波程序clear;clc;Fs=1000;nfft=2048;t1=0:1/Fs:2.047;A=normrnd(0,1,1,2048);N=wgn(1,2048,2);f1=390;f2=410;wc1=2*f1/Fs;wc2=2*f2/Fs;wc2=2*f2/sf;%歸一化頻率f0f0=[0wc1-0.05wc1wc2wc2+0.051];B=[001100];%設立帶通和帶阻weigh=[111];%設立帶通和帶阻權重b=remez(50,f0,B,weigh);%傳函分子D=filter(b,1,N);x=sin(2*pi*t1*100)+1.5*sin(2*pi*t1*300)+A.*sin(2*pi*t1*200)+D+N;x1=sin(2*pi*t1*100)+1.5*sin(2*pi*t1*300)+A.*sin(2*pi*t1*200);a1=-1.352;a2=1.338;a3=-0.662;a4=0.240;A=[-a1-a2-a3-a4;1000;0100;0010];%狀態(tài)轉移矩陣H=[1000];%觀察矩陣Q=[1000;0000;0000;0000];%狀態(tài)噪聲方差R=1;%觀察噪聲方差陣X(:,1)=[x(4);x(3);x(2);x(1)];p(:,:,1)=[10000;0100;0010;0001];%一步預測誤差方針%開始濾波fork=2:nfftp1(:,:,k)=A*p(:,:,k-1)*A'+Q;%p1(:,:,k)即是一步預測誤差的自有關矩陣，它是4*4的矩陣，取不同的k值就構成了一種三維矩陣K(:,k)=p1(:,:,k)*H'/(H*p1(:,:,k)*H'+R);%K(:,:,k)是增益矩陣,對于固定的k值它是4*1矩陣，取不同的k值就是三維矩陣X(:,k)=A*X(:,k-1)+K(:,k)*[x(k)-H*A*X(:,k-1)];%X(:,k)是預計值,4*1矩陣p(:,:,k)=p1(:,:,k)-K(:,k)*H*p1(:,:,k);%p(:,:,k)是預計誤差的自有關矩陣，4*4矩陣的三維矩陣end%結束一次濾波%繪圖t=1:nfft;figure(2);plot(t,x1,'k-',t,x,'r-',t,X(1,:),'b-.');title('卡曼濾波去噪')legend('真實軌跡','觀察樣本','預計軌跡');gridon;圖5Kalman濾波去噪圖(4)運用Wigner-Ville分布分析信號的時頻特性MATLAB程序clear;clc;Fs=1000;nfft=2049;t1=0:1/Fs:2.048;A=normrnd(0,1,1,2049);N=wgn(1,2049,2);f1=390;f2=410;wc1=2*f1/Fs;wc2=2*f2/Fs;%歸一化頻率f0f0=[0wc1-0.05wc1wc2wc2+0.051];B=[001100]