4.信號(hào)的函數(shù)體現(xiàn)式為：，其中，為一隨時(shí)間變化的隨機(jī)過(guò)程，為通過(guò)390—410Hz帶通濾波器后的高斯白噪聲，為高斯白噪聲，采樣頻率為1kHz，采樣時(shí)間為2.048s。（1）運(yùn)用當(dāng)代信號(hào)解決的知識(shí)進(jìn)行信號(hào)譜預(yù)計(jì)；（2）運(yùn)用當(dāng)代信號(hào)解決知識(shí)進(jìn)行信號(hào)的頻率提??；（3）分別運(yùn)用Winner濾波和Kalman濾波進(jìn)行去噪；（4）運(yùn)用Wigner-Ville分布分析信號(hào)的時(shí)頻特性。（1）：運(yùn)用當(dāng)代信號(hào)解決的知識(shí)進(jìn)行信號(hào)譜預(yù)計(jì)：典型譜預(yù)計(jì)中兩種重要的辦法為直接法和間接法，其中間接法則先根據(jù)N個(gè)樣本數(shù)據(jù)的樣本自有關(guān)函數(shù) （4.1）其中，且。計(jì)算樣本自有關(guān)函數(shù)的Fourier變換，得到功率譜 （4.2）周期圖辦法預(yù)計(jì)的功率譜為有偏預(yù)計(jì)，可通過(guò)加窗來(lái)減少其偏差。定義為 （4.3）式中 （4.4）式中，是窗函數(shù)的Fourier變換。功率譜預(yù)計(jì)程序?yàn)椋篶learclccloseallhiddensf=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];%設(shè)立帶通和帶阻weigh=[111];%設(shè)立帶通和帶阻權(quán)重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功率譜預(yù)計(jì)成果圖(2).信號(hào)頻率的提取用離散傅立葉算法離散傅立葉算法程序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];%設(shè)立帶通和帶阻weigh=[111];%設(shè)立帶通和帶阻權(quán)重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('原始信號(hào)');axis([02.047-68]);plot(f,y1);title('fft頻率提取');axis([050001000]);xlabel('f/Hz');gridon;圖4.2原始信號(hào)時(shí)域圖圖4.3信號(hào)頻譜(3)分別運(yùn)用Winner濾波和Kalman濾波進(jìn)行去噪；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];%設(shè)立帶通和帶阻weigh=[111];%設(shè)立帶通和帶阻權(quán)重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';%獲得抱負(fù)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%繪圖比較預(yù)計(jì)濾波器與實(shí)際濾波器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)');%比較抱負(fù)輸出與實(shí)際輸出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];%設(shè)立帶通和帶阻weigh=[111];%設(shè)立帶通和帶阻權(quán)重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)轉(zhuǎn)移矩陣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];%一步預(yù)測(cè)誤差方針%開始濾波fork=2:nfftp1(:,:,k)=A*p(:,:,k-1)*A'+Q;%p1(:,:,k)即是一步預(yù)測(cè)誤差的自有關(guān)矩陣，它是4*4的矩陣，取不同的k值就構(gòu)成了一種三維矩陣K(:,k)=p1(:,:,k)*H'/(H*p1(:,:,k)*H'+R);%K(:,:,k)是增益矩陣,對(duì)于固定的k值它是4*1矩陣，取不同的k值就是三維矩陣X(:,k)=A*X(:,k-1)+K(:,k)*[x(k)-H*A*X(:,k-1)];%X(:,k)是預(yù)計(jì)值,4*1矩陣p(:,:,k)=p1(:,:,k)-K(:,k)*H*p1(:,:,k);%p(:,:,k)是預(yù)計(jì)誤差的自有關(guān)矩陣，4*4矩陣的三維矩陣end%結(jié)束一次濾波%繪圖t=1:nfft;figure(2);plot(t,x1,'k-',t,x,'r-',t,X(1,:),'b-.');title('卡曼濾波去噪')legend('真實(shí)軌跡','觀察樣本','預(yù)計(jì)軌跡');gridon;圖5Kalman濾波去噪圖(4)運(yùn)用Wigner-Ville分布分析信號(hào)的時(shí)頻特性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]