wiener滤波程序

1、基于最小均方误差(MMSE)估计的因果维纳滤波的实现一功能简介 基于最小均方误差(MMSE)估计的因果维纳滤波的Matlab实现，用莱文森-德宾(Levinson-Durbin)算法求解维纳-霍夫方程(Yule-wa1ker)方程，得到滤波器系数，进行维纳滤波。二维纳滤波简介信号处理的实际问题，常常是要解决在噪声中提取信号的问题，因此，我们需要寻找一种所谓有最佳线性过滤特性的滤波器，这种滤波器当信号与噪声同时输入时，在输出端能将信号尽可能精确地重现出来，而噪声却受到最大抑制。维纳(Wiener)滤波就是用来解决这样一类从噪声中提取信号问题的一种过滤(或滤波)方法。一个线性系统，如果它的单位样本

2、响应为h(n)，当输入一个随机信号x(n)，且其中s(n)表示信号，表示噪声，则输出y(n)为我们希望x(n)通过线性系统h(n)后得到的y(n)尽量接近于s(n)，因此称y(n)为s(n)的估计值，用表示，即维纳滤波器的输入输出关系如上图所示。这个线性系统称为对于的一种估计器。如果我们以分别表示信号的真值与估计值，而用e(n)表示它们之间的误差，即显然，e(n)可能是正的，也可能是负的，并且它是一个随机变量。因此，用它的均方值来表达误差是合理的，所谓均方误差最小即它的平方的统计平均值最小：最小已知希望输出为：误差为：均方误差为：上式对求导得到：进一步得：从而有：于是就得到N个线性方程：写成矩

3、阵形式为：简化形式：其中：是滤波器的系数是互相关序列是自相关矩阵由上可见，设计维纳滤波器的过程就是寻求在最小均方误差下滤波器的单位脉冲响应或传递函数的表达式，其实质就是解维纳霍夫（WienerHopf）方程。另外，设计维纳滤波器要求已知信号与噪声的相关函数。三程序求解过程由上述可见，本程序实现的关键是在已知输入信号的自相关函数和输入信号和理想输出信号的互相关函数的情况下，求解维纳霍夫（WienerHopf）方程，从而得到滤波器系数，再进行维纳滤波。求解步骤：1. 初始化值2. 对于，进行如下计算：3滤波器系数为：4利用上面的得到的滤波器对输入信号进行维纳滤波，得到输出信号。四函数说明函数使用方

4、法：y=wienerfilter(x,Rxx,Rxd,M)参数说明：x是输入信号，Rxx是输入信号的自相关向量，Rxd是输入信号和理想信号的的互相关向量，M是维纳滤波器的长度，输出y是输入信号通过维纳滤波器进行维纳滤波后的输出。具体程序见Matlab的.m文件。五程序示例加载Matlab中的语音数据handel，人为地加入高斯白噪声，分别计算加入噪声后信号的自相关和加入噪声后信号和理想信号的互相关，取滤波器的长度为M=500，将以上参数代入函数中进行维纳滤波，得到输出。程序如下：load handel %加载语音信号d=y; d=d*8; %增强语音信号强度d=d'fq=fft(d,8

5、192); %进行傅立叶变换得到语音信号频频subplot(3,1,1);f=Fs*(0:4095)/8192;plot(f,abs(fq(1:4096); %画出频谱图title('原始语音信号的频域图形');xlabel('频率 f');ylabel('FFT');m,n=size(d);x_noise=randn(1,n); %（0，1）分布的高斯白噪声x=d+x_noise; %加入噪声后的语音信号fq=fft(x,8192); %对加入噪声后的信号进行傅立叶变换，看其频谱变化subplot(3,1,2);plot(f,abs(fq(1:

6、4096); %画出加入噪声后信号的频谱图title('加入噪声后语音信号的频域图形');xlabel('频率 f');ylabel('FFT');yyhxcorr=xcorr(x(1:4096); %求取信号的信号的自相关函数size(yyhxcorr); A=yyhxcorr(4096:4595);yyhdcorr=xcorr(d(1:4096),x(1:4096); %求取信号和理想信号的互相关函数size(yyhdcorr);B=yyhdcorr(4096:4595);M=500;yyhresult=wienerfilter(x,A,B,

7、M); %进行维纳滤波yyhresult=yyhresult(300:8192+299);fq=fft(yyhresult); %对维纳滤波的结果进行傅立叶变换，看其频谱变化subplot(3,1,3); f=Fs*(0:4095)/8192;plot(f,abs(fq(1:4096); %画出维纳滤波后信号的频谱图title('经过维纳滤波后语音信号的频域图形');xlabel('频率 f');ylabel('FFT');求出的频谱图如下所示：由上述结果可见，经过维纳滤波后信号的噪声减弱，信噪比提高。yyhwiener.mload handel

8、 %sound(5*y,Fs) %原始语音信号,Fs=9192 d=y; d=d*8; d=d' fq=fft(d,8192); subplot(3,1,1); f=Fs*(0:4095)/8192; plot(f,abs(fq(1:4096); title('原始语音信号的频域图形'); xlabel('频率 f'); ylabel('FFT'); m,n=size(d); x_noise=randn(1,n); x=d+x_noise; %加入噪声后的语音信号，噪声为(0,1)分布的高斯白噪声 fq=fft(x,8192); subp

9、lot(3,1,2); plot(f,abs(fq(1:4096); title('加入噪声后语音信号的频域图形'); xlabel('频率 f'); ylabel('FFT'); yyhxcorr=xcorr(x(1:4096); size(yyhxcorr); A=yyhxcorr(4096:4595); yyhdcorr=xcorr(d(1:4096),x(1:4096); size(yyhdcorr); B=yyhdcorr(4096:4595); M=500; yyhresult=wienerfilter(x,A,B,M); %用维纳滤

10、波进行去噪 yyhresult=yyhresult(300:8192+299); fq=fft(yyhresult); subplot(3,1,3); f=Fs*(0:4095)/8192; plot(f,abs(fq(1:4096); title('经过维纳滤波后语音信号的频域图形'); xlabel('频率 f'); ylabel('FFT');wienerfilter.mfunction y=wienerfilter(x,Rxx,Rxd,N) %进行维纳滤波 %x是输入信号，Rxx是输入信号的自相关向量 %Rxx是输入信号和理想信号的的互相

11、关向量，N是维纳滤波器的长度 %输出y是输入信号通过维纳滤波器进行维纳滤波后的输出 h=yulewalker(Rxx,Rxd,N);%求解维纳滤波器系数 t=conv(x,h); %进行滤波 Lh=length(h); %得到滤波器的长度 Lx=length(x); %得到输入信号的长度 y=t(double(uint16(Lh/2):Lx+double(uint16(Lh/2)-1);%输出序列y的长度和输入序列x的长度相同 %以下是维纳滤波器系数的求解 function h=yulewalker(A,B,M) %求解Yule-Walker方程 %A是接收信号的自相关向量为 Rxx(0),R

12、xx(1),.,Rxx(M-1) %B是接收信号和没有噪声干扰信号的互相关向量为 Rxd(0),Rxd(1),.,Rxd(M-1) %M是滤波器的长度 %h保存滤波器的系数 %例如A=6,5,4,3,2,1;B=100,90,120,50,80,200;M=6; %求解出h=26.4286 -20.0000 50.0000 -50.0000 -45.0000 81.4286' T1=zeros(1,M);%T1存放中间方程的解向量 T2=zeros(1,M);%T2存放中间方程的解向量 T1(1)=B(1)/A(1); T2(1)=A(2)/A(1); X=zeros(1,M); fo

13、r i=2:M-1 temp1=0; temp2=0; for j=1:i-1 temp1=temp1+A(i-j+1)*T1(j); temp2=temp2+A(i-j+1)*T2(j); end X(i)=(B(i)-temp1)/(A(1)-temp2); for j=1:i-1 X(j)=T1(j)-X(i)*T2(j); end for j=1:i T1(j)=X(j); end temp1=0; temp2=0; for j=1:i-1 temp1=temp1+A(j+1)*T2(j); temp2=temp2+A(j+1)*T2(i-j); end X(1)=(A(i+1)-temp1)/(A(1)-temp2); for j=2:i X(j)=T2(j-1)-X(1)*T2(