白化滤波器-matlab-程序

白化滤波器原理在统计信号处理中，往往会遇到等待处理的随机信号是非白色的，例如云雨、海浪、地物反射的杂乱回波等，它们的功率谱即使在信号通带内也非均匀分布。这样会给问题的解决带来困难。克服这一困难的措施之一是对色噪声进行白化处理。主要内容是设计一个稳定的线性滤波器，将输入的有色噪声变成输出的白噪声。在这里，我们就对一般的具有功率谱)(xG的平稳随机过程X(t)白化处理问题进行讨论。为了具体的进行分析和计算，假设)(xG可以表达成有理数的形式，即))......(())......(()(112mnxZZaGmnZ其中分子、分母为多项式。这个假设对于通常见到的功率谱是很近似的，而且有可行的方法用有理数去逼近任意的功率谱密度。由于)(xG是功率谱，它的平稳随机过程相关函数的傅里叶变换具有非负的实函数和偶函数的性质。这些性质必然在其有理函数的表示式中体现出来，特别是，)(xG的零、极点的分布和数量会具有若干个特点。由于)(xG是实函数，因此有：)()(*xxGG，2a是实数，)(xG的零、极点是共轭成对的。从而也可以把)(xG的表示式写成如下形式：))......(())......(())......(())......(()(1111lklkxjjjjajjjjaG把开拓到复平面s中去，另js。用s代替j就可以把函数)(xG扩大到整个复平面。)(xG的零、极点必将对称于轴，如图13所示：图13由于)(xG是偶函数，因此不难判断，)(xG的零、极点是象限对称的，从而对于j轴也是对随机信号分析实验-2-称的。由于0)(xG，因此分子的虚根必然是偶数，否则)(xG会出现负值。这就是说j轴上的零、极点必将成对的出现。由于)(xG是可积的，因此分子的阶数不能大于分母的阶数，这就是说零点总数不会大于极点总数，而且分母不可能有虚根，这意味着j轴上没有极点。综合上述情况，在s平面的零、极点的可能位置如上图所示：令：))......(())......(()(11lkxjjjjaG))......(())......(()(11lkxjjjjaG则有)()()()()(*xxxxxGGGGG其中)(xG代表零、极点均在s左平面的部分，)(xG代表零、极点均在s右平面的部分。若在j轴上有零点的话，必是成对的。则将一个放在)(xG内，将另一个放在)(xG内。实质上，)(xG对应的时域函数在负时间域为零，而)(xG对应的时域函数在正时间域为零。根据上述的讨论，可以求得白化滤波器)(1H的解析式为：)(1|)(|21xGH由于)()()()(|)(|11*1121HHHHH)()()()()(xxxxxGGGGG故得：)(1)(1xGH若运用傅里叶变换进行分析计算，以s代替j，可得白化滤波器公式：)(1)(1sGsHx其中js我们知道，)(1H的傅里叶反变换是白化滤波器在时域的单位冲击响应)(1th，)(xG零、极点随机信号分析实验-3-在s左半平面，因此)(1xG的零、极点也是在s左半平面。故它对应的时域函数)(1th在负时域时为零，也就是说，上述白化滤波器是物理可实现的。白化滤波器的设计方法是：首先计算色噪声自相关函数，根据色噪声的自相关函数，计算出色噪声的功率谱（色噪声的自相关函数和功率谱构成一对傅里叶变换对），然后根据公式)(1|)(|21xGH程序%*******色噪声的产生*******************%Fs=44100;[x1,Fs]=wavread('E:\matlab\work\混合信号\色噪声.wav');L1=length(x1)/10;x=x1(1:44100);l1=0:L1-1;t=l1/Fs;figure(1);plot(t,x,'-r');title('色噪声');%******统计色噪声*******************E=mean(x);%色噪声的均值E=-0.0054S=var(x);%色噪声的方差S=0.0324%**************求色噪声概率密度函数*****************%eachi=linspace(min(x),max(x),42);yyi=hist(x,eachi);%计算各个区间的个数yyi=yyi/length(x);%对各个区间的个数归一化处理figure(8);%绘制色噪声的概率密度函数plot(eachi,yyi,'-k')title('色噪声的概率密度函数')%***************色噪声自相关函数*****%Rx=xcorr(x,x);%色噪声的自相关函数Rxtau=(-L1+1:L1-1)/Fs;figure(2);plot(tau,Rx,'-r')title('色噪声的自相关函数');%色噪声的自相关函数波形xlabel('\tau'),ylabel('R_x(\tau)');随机信号分析实验-4-gridon;holdon;%***************色噪声功率谱密度*****%R=fft(Rx);%自相关函数的傅里叶变换即是功率谱密度cm=abs(R);fl=(0:length(R)-1)*44100/length(R);figure(3)plot(fl(1:length(fl)/2),cm(1:length(fl)/2),'-b')title('色噪声的功率谱')holdon;gridon;Yz=length(z);Lz=0:Yz-1;Tz=Lz/Fs;figure(9)plot(Tz,z,’-r’);title(‘白噪声’)%******白化滤波器的产生****************%k=sqrt(cm);c=1./k;figure(4);plot(fl(1:length(fl)/2),c(1:length(fl)/2),'-b')title('白化滤波器滤波特性')%****************统计白噪声**************%ci=ifft(c);%白化滤波器在时域特性h（t）z=conv(x,ci);%（色噪声与h（t）的卷积，输出白噪声）Eo=mean(z);%白噪声的均值Eo=-7.5587e-006So=var(z);%白噪声的方差So=6.2662e-005%**************求白噪声概率密度函数*****************%eacho=linspace(min(z),max(z),42);yyo=hist(z,eacho);%计算各个区间的个数yyo=yyo/length(z);%对各个区间的个数归一化处理figure(7);%绘制白噪声的概率密度函数plot(eacho,yyo,'-k')title('白噪声的概率密度函数')随机信号分析实验-5-%***************输出信号自相关函数*****%Rz=xcorr(z,z);%输出白化噪声的自相关函数Rztau2=(-length(z)+1:length(z)-1)/Fs;figure(5);plot(tau2,Rz,'-r')title('输出白化噪声的自相关函数');%输出白化信号的自相关函数波形xlabel('\tau2'),ylabel('R_z(\tau2)');gridon;holdon;%***************输出信号功率谱密度*****%Ro=fft(Rz);%自相关函数的傅里叶变换即是功率谱密度cmo=abs(Ro);f2=(0:length(Ro)-1)*44100/length(Ro);figure(6);plot(f2(1:length(f2)/2),cmo(1:length(f2)/2),'-b');title('输出信号的功率谱');holdon;gridon;