现代信号处置

4.信号的函数体现式为：，其中，为一随时间变化的随机过程，为通过390—410Hz带通滤波器后的高斯白噪声，为高斯白噪声，采样频率为1kHz，采样时间为2.048s。（1）运用当代信号解决的知识进行信号谱预计；（2）运用当代信号解决知识进行信号的频率提取；（3）分别运用Winner滤波和Kalman滤波进行去噪；（4）运用Wigner-Ville分布分析信号的时频特性。（1）：运用当代信号解决的知识进行信号谱预计：典型谱预计中两种重要的办法为直接法和间接法，其中间接法则先根据N个样本数据的样本自有关函数 （4.1）其中，且。计算样本自有关函数的Fourier变换，得到功率谱 （4.2）周期图办法预计的功率谱为有偏预计，可通过加窗来减少其偏差。定义为 （4.3）式中 （4.4）式中，是窗函数的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;%维纳滤波器阶数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滤波器系数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];%状态转移矩阵H=[1000];%观察矩阵Q=[1000;0000;0000;0000];%状态噪声方差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]