维纳滤波器实验

维纳滤波器实验 完整的维纳滤波器Matlab源程序 tp://www.voanews.com/english/news/ clear;clc; %输入信号 A=1; %信号的幅值 f=1000; %信号的频率 fs=10^5; %采样频率 t=(0:999); %采样点 Mlag=100; %相关函数长度变量 x=A*(sin(2*pi*f*t/fs)+sin(4*pi*f*t/fs)+sin(6*pi*f*t/fs)); %输入正弦波信号 xmean=mean(x); %正弦波信号均值 xvar=var(x,1); %正弦波信号方差 xn=awgn(x,5); %给正弦波信号加入信噪比为20dB的高斯白噪声 figure(1) plot(t,xn) %绘制输入信号图像 title('输入信号图像') xlabel('x轴单位:t/s','color','b') ylabel('y轴单位:f/HZ','color','b') xnmean=mean(xn) %计算输入信号均值 xnms=mean(xn.^2) %计算输入信号均方值 xnvar=var(xn,1) %计算输入信号方差 Rxn=xcorr(xn,Mlag,'biased'); %计算输入信号自相关函数 figure(2) subplot(221) plot((-Mlag:Mlag),Rxn) %绘制自相关函数图像 title('输入信号自相关函数图像') [f,xi]=ksdensity(xn); %计算输入信号的概率密度，f为样本点xi处的概率密度 subplot(222) plot(xi,f) %绘制概率密度图像 title('输入信号概率密度图像') X=fft(xn); %计算输入信号序列的快速离散傅里叶变换 Px=X.*conj(X)/600; %计算信号频谱 subplot(223) semilogy(t,Px) %绘制在半对数坐标系下频谱图像 title('输入信号在半对数坐标系下频谱图像') xlabel('x轴单位:w/rad','color','b') ylabel('y轴单位:w/HZ','color','b') pxx=periodogram(xn); %计算输入信号的功率谱密度 subplot(224) semilogy(pxx) %绘制在半对数坐标系下功率谱密度图像 title('输入信号在半对数坐标系下功率谱密度图像') xlabel('x轴单位:w/rad','color','b') ylabel('y轴单位:w/HZ','color','b') %fir滤波 wp=3000; %通带截止频率 ws=4000; %阻带截止频率 DB=ws-wp; %过渡带宽度 N0=ceil(6.6*pi/DB); M=N0+mod(N0+1,2); %计算fir滤波器阶数 wc=(wp+ws)/2/pi; %计算理想低通滤波器通带截止频率(关于π归一化) hn=fir1(M,wc); %调用fir1计算FIRDF的h(n) y1n=filter(hn,1,xn); %将输入信号通过fir滤波器 figure(3) plot(y1n) %绘制经过fir滤波器后信号图像 title('经过fir滤波器后信号图像') xlabel('x轴单位:f/HZ','color','b') ylabel('y轴单位:A/V','color','b') y1nmean=mean(y1n) %计算经过fir滤波器后信号均值 y1nms=mean(y1n.^2) %计算经过fir滤波器后信号均方值 y1nvar=var(y1n,1) %计算经过fir滤波器后信号方差 Ry1n=xcorr(y1n,Mlag,'biased'); %计算经过fir滤波器后信号自相关函数 figure(4) subplot(221) plot((-Mlag:Mlag),Ry1n) %绘制自相关函数图像 title('经过fir滤波器后信号自相关函数图像') [f,y1i]=ksdensity(y1n); %计算经过fir滤波器后信 号的概率密度，f为样本点xi处的概率密度 subplot(222) plot(y1i,f) %绘制概率密度图像 title('经过fir滤波器后信号概率密度图像') Y1=fft(y1n); %计算经过fir滤波器后信号序列的快速离散傅里叶变换 Py1=Y1.*conj(Y1)/600; %计算信号频谱 subplot(223) semilogy(t,Py1) %绘制在半对数坐标系下频谱图像 title('经过fir滤波器后信号在半对数坐标系下频谱图像') xlabel('x轴单位:w/rad','color','b') ylabel('y轴单位:w/HZ','color','b') py1n=periodogram(y1n); %计算经过fir滤波器后信号的功率谱密度 subplot(224) semilogy(py1n) %绘制在半对数坐标系下功率谱密度图像 title('经过fir滤波器后信号在半对数坐标系下功率谱密度图像') xlabel('x轴单位:w/rad','color','b') ylabel('y轴单位:w/HZ','color','b') %维纳滤波 N=100; %维纳滤波器长度 Rxnx=xcorr(xn,x,Mlag,'biased'); %产生输入信号与原始信号的互相关函数 rxnx=zeros(N,1); rxnx(:)=Rxnx(101:101+N-1); Rxx=zeros(N,N); %产生输入信号自相关矩阵 Rxx=diag(Rxn(101)*ones(1,N)); for i=2:N c=Rxn(101+i)*ones(1,N+1-i); Rxx=Rxx+diag(c,i-1)+diag(c,-i+1); end Rxx; h=zeros(N,1); h=inv(Rxx)*rxnx; %计算维纳滤波器的h(n) yn=filter(h,1,xn); %将输入信号通过维纳滤波器 figure(5) plot(yn) %绘制经过维纳滤波器后信号图像 title('经过维纳滤波器后信号信号图像') xlabel('x轴单位:f/HZ','color','b') ylabel('y轴单位:A/V','color','b') ynmean=mean(yn) %计算经过维纳滤波器后信号均值 ynms=mean(yn.^2) %计算经过维纳滤波器后信号均方值 ynvar=var(yn,1) %计算经过维纳滤波器后信号方差 Ryn=xcorr(yn,Mlag,'biased'); %计算经过维纳滤波器后信号自相关函数 figure(6) subplot(221) plot((-Mlag:Mlag),Ryn) %绘制自相关函数图像 title('经过维纳滤波器后信号自相关函数图像') [f,yi]=ksdensity(yn); %计算经过维纳滤波器后信号的概率密度，f为样本点xi处的概率密度 subplot(222) plot(yi,f) %绘制概率密度图像 title('经过维纳滤波器后信号概率密度图像') Y=fft(yn); %计算经过维纳滤波器后信号序列的快速离散傅里叶变换 Py=Y.*conj(Y)/600; %计算信号频谱 subplot(223) semilogy(t,Py) %绘制在半对数坐标系下频谱图像 title('经过维纳滤波器后信号在半对数坐标系下频谱图像') xlabel('x轴单位:w/rad','color','b') ylabel('y轴单位:w/HZ','color','b') pyn=periodogram(yn); %计算经过维纳滤波器后信号的功率谱密度 subplot(224) semilogy(pyn) %绘制在半对数坐标系下功率谱密度图像 title('经过维纳滤波器后信号在半对数坐标系下功率谱密度图像') xlabel('x轴单位:w/rad','color','b') ylabel('y轴单位:w/HZ','color','b')