白噪声中正弦组合的功率谱分析.doc

通信工程 孙大江 郑州大学 白噪声中正弦组合的功率谱分析 专业： 通信工程(1)学生姓名： 孙大江学号: 20092420128指导教师： 陈恩庆 完成时间： 2020年1月29日 摘要 本文分别选取古典谱估计,现代参数方程谱估计,特征分解法谱估计中的各一种对淹没在白噪声中的正弦信号进行功率谱分析,并进行了相应的matlab仿真,仿真程序和结果见附录. 关键词:功率谱,自相关,估计 Abstract This paper selected the classical spectrum estimation, parameter equation of modern spectral estimation, feature decomposition spectral estimation to analysis the sine signal concealed in the white noise , and conducted the implementation of the corresponding matlab simulation ,the Simulation program and results are Save on the appendix. Key words: power spectrum, autocorrelation, estimation 一.课题简介 由于白噪声中的正弦组合是最常见的随机过程,而估计淹没在噪声中的正弦波频率即是信号处理中最有实际应用价值的技术之一,也是测试所有谱估计性能的基础.所以在本课题中以信号 为例进行研究,并分别以经典谱估计中的自相关法(间接法),和现代谱估计AR模型的Yule-Walker法(自相关法),以及分离信号子空间和噪声子空间的特征分解法(MUSIC)分别进行功率谱估计.最后再根据不同方法得到的不同结果进行分析.二.信号分析的演变在进行课题研究之前,我们首先明确一下为什么要进行信号谱分析.最初的信号研究知识局限在时域,但是随着科技的进步,时域分析已经完全满足不了科研需要,于是人们开始把目光转向了频域. 最初,根据任何信号f(t)都可以在特定区间内展开成由完备正交集构成的无穷级数的理论,我们可以吧一个周期信号分解为有不同频率谐波组成的形式,并且分别在三角集和指数集中得到三角型傅里叶级数和指数型傅里叶级数,由于指数型运算的简便性,人们选择了指数型运算进行了深入发展,即以为振幅,以为相角分别得到周期信号的频域幅度谱和相位谱对于非周期信号,由于相当于周期T区域无穷,所以基波角频率区域无穷小,个分量的幅度也相应的趋于无穷小,不能在进行幅度谱分析了,故引出新的概念频谱密度 即单位频率的振幅,也就是我们通常所说的频谱 ,这里引出了傅里叶变换,即根据傅里叶变换,我们可以进行频谱分析了 下一步,引入能量谱的概念,信号能量,其中为单位频率的信号能量,也称为能量密度函数,即我们所说的能量谱 , 第四步,与能量谱类似引入功率谱概念 , 由于随机信号不能用频谱表示(能量无限),但是由维纳辛钦定理,功率有限信号的功率谱函数语气自相关函数是一对傅里叶变换,所以我们可以由信号的自相关函数求得其功率谱,然后利用功率谱来描述随机信号的频域特性.三.本文中用到的三种谱分析方法.古典谱估计之相关法 相关法谱估计是以相关函数为媒介来计算功率谱,又叫做间接法它的理论基础是维纳-辛钦定理,其具体步骤如下: 第一步,由获得的N点数据构成的有限长序列xn(n)来估计自相关函数,即:第二步,由自相关函数的傅里叶变换求功率谱,即以上两步经历了两次截断,一次是估计 时仅利用了x(n)的N个观测值,这相当于对x(n)加矩形窗截断.该窗是加在数据上的,一般称为加数据窗.另一次是估计时仅利用了从-(M-1)到(M-1)的,这相当于对加矩形窗截断,将截成(2M-1)长,这称为加延时窗.式中和分别表示对和的估值.由于MN,使得上式的运算量不是很大,因此在FFT问世之前,相关法是最常用的谱估计方法. FFT问世以后,情况有所变化,这里只进行基本思想介绍,不再进行matlab程序仿真. 因为阶段后的可视作能量信号,由相关卷积定理对上式两边取(2N-1)点DFT,并将时域卷积变为频域乘积,有于是可得用FFT求方案如下:(1) 对N长充(N-1)个零,成为(2N-1)长的.(2) 求(2N-1)点的FFT,得:(3) 求(4) 求(2N-1)点的IFFT,得: 上面的相关运算中,充零的目的是为了能用圆周卷积代替线性卷积,以便采用快速卷计算法.,参数模型法谱估计之levinson-durbin快速递推法 Levinson-durbin递推算法是解尤勒-沃克方程的快速有效的算法,这种方法利用自相关矩阵具有的一些好的性质,是运算量大大减少,在介绍Levinson-durbin递推算法之前,先大概了解一下尤勒-沃克方程: 如前所述,P阶模型的系统函数为 可以看出，阶模型有个待定系数，由上式，可得白噪声激励得到的系统输出 可以理解为，用时刻之前的个值的线性组合来预测时刻的值预测误差为在均方误差最小准则下，组合系数的选择应使预测误差的均方值最小经过一系列的运算，最终可以得到模型的正则方程也就是尤勒沃克（）方程在上面方程中令，得一阶模型的尤勒-沃克方程为可以解得进而求得 在一阶的基础上进行递推，将阶次为时的第个系数定义为反射系数，用表示，于是可以将计算阶模型参数的Levinson-durbin递推算法表示如下：式中，.特征分解法谱估计之MUSIC方法前面分别讨论了古典法和参数模型法的一个例子，下面讨论用特征分解法对白噪声中的正弦波频率进行谱估计特征分解技术的主要思想是，把自相关矩阵中的信息空间分成两个子空间，即信号子空间和噪声子空间，这两个子空间中的矢量的函数在正弦波频率上有尖锐的峰，据此可以估计正弦波的频率输入信号为 其向量形式为 由向量可以求出自相关矩阵为 其中，分别为信号自相关矩阵和噪声自相关矩阵也就是说，数据自相关矩阵可以分解为信号自相关矩阵和噪声自相关矩阵之和和都是N行N列方阵,秩分别为M和N,即 .如果已知或者其估计,则可以通过分解由得到,再进一步分解由得到,从而求出正弦波频率的估计.省略分解的过程,定义N个特征向量中非零的为主特征向量为噪声特征向量,因为信号子空间与噪声子空间相互正交,所以信号向量与噪声子空间的向量都正交,与它们的线性组合也正交,有 式中为M个正弦信号的频率.令 当时,应有 即 因此可以定义一中类似于功率谱的函数若令,则有上式取峰值的M个就是M个正弦信号频率的估计.理论上这M个应使上式所示的函数值为无限大,但由于存在估计误差,所以为有限值,但呈现尖锐的峰值.也就是说,的M个最大值所对应的频率就是正弦信号频率的估计.4. 程序和结果时域图如下:(1)间接法程序clear;n=1:1000;x=sqrt(20)*sin(2*pi*0.2*n)+sqrt(2)*sin(2*pi*0.3*n)+randn(size(n);N=1000;for h=1:N R(h)=0;endfor m=1:N s=0; for i=1:N-(m-1) s=s+x(i)*x(i+m-1); end R(m)=s/N; %估计自相关函数,P89页4.1.1end w=0:0.001:0.5; Sx=0;M=1000;for k=1:(M-1) Sx=Sx+R(k)*exp(-j*w*k)+R(k)*exp(j*w*k);%由自相关函数的傅里叶变换求 功率谱,P89页4.1.2endSx=log10(abs(Sx); plot(w,abs(Sx);title(间接法谱估计);grid on;(2)Levinson法程序clc;clear allfs=100; %设采样频率为100，N=150 ; %数据长度,改变数据长度会导致分辨率的变化；n=0:0.01:2; %变化范围x=sqrt(20)*sin(2*pi*0.2*n)+sqrt(10)*sin(2*pi*0.3*n)+randn(size(n);for m=1:N %接下来的循环计算自相关函数的无偏估计,仍使用P89也方法 R(m)=0;endfor m=1:N s=0; for n=1:N-(m-1) s=s+x(n)*x(n+m-1); end R(m)=s/N;end for M=1:N-5 %接下来的循环定义初值-阶数M要小于数据长度N a(M,M)=0; FPE(M)=0; P(M)=0;end %接下来初始化：由Levinson公式a(1,1)=-R(2)/R(1); %计算一阶模型参数a(1,1),见P66页公式3.3.15；P(1)=R(1)+a(1,1)*R(2); %计算一阶模型的误差功率P(1),p66页底部,误差功率P=G2；FPE(1)=P(1)*(N+1+1)/(N-1-1) %计算一阶模型的最终预测误差准则函数FPE(1)；sum=0; %接下来进行递推for M=2:N-5 for k=1:M-1 sum=sum+a(M-1,k)*R(M-k+1); end a(M,M)=-(R(M+1)+sum)/P(M-1); %P67页计算反射系数公式(3.3.16a) for k=1:M-1 a(M,k)=a(M-1,k)+a(M,M)*a(M-1,M-k); %由反射系数通过levinson算法递推求系统函数 系数a(3.3.16b) end P(M)=(1-(abs(a(k,k)2)*P(M-1); %由反射系数求预测误差功率,即G2(3.3.16c) FPE(M)=P(M)*(N+M+1)/(N-M-1); %计算模型的最终预测误差准则函数FPE(M)； sum=0; %sum归零end %确定阶数:求出使FPE最小时的阶数M1；min=FPE(1); %确定比较对象为FPE(1)for M=2:N-5 %依次比较 if FPE(M)min min=FPE(M); M1=M; endend disp(阶数M为：) %确定阶数disp(M1);disp(模型参数a为：); %由得到的阶数只截取参数a从1到M1 for i=1:M1 disp(a(M1,i);enddisp(误差功率为：); %FPE最小时误差功率最小,得到需要的误差功率,最终得到系统函数disp(P(M1); %AR模型参数确定后计算出功率谱Z=0;W=0:0.01:pi; %功率谱以2*pi为周期，又信号为实信号，只需输出0到PI即可；for k=1:M1 Z=Z+(a(M1,k).*exp(-j*k*W); %计算AR模型分母 endPSD=P(M1)./(abs(1+Z).2); %计算功率谱,系统函数H(z)平方即为模型能量谱在数值上 %等于功率谱,(因为白噪声激励,自功率谱为1)P61(3.2.8)F=W*fs/(pi*200); %角频率坐标换算成频率坐标plot(F,abs(PSD); title(Levison Durbin谱估计);grid on(3)特征分解法谱估计clear alln=1:20; %信号长度N=20; %自相关矩阵的秩x(n)=sqrt(20)*sin(2*pi*0.2*n)+sqrt(2)*sin(2*pi*0.3*n)+randn(size(n);X=zeros(N,1);for i=1:N X(i,1)=x(i); %对全零矩阵各元素重新赋值生成X信号矩阵 N行1 列endrxx=zeros(N,N); %定义一个全零矩阵 N行N列for i=1:N; rxx=rxx+X(i,:)*(X(i,:); %矩阵X乘以X的共轭转置,得到N行N列矩阵endrxx=rxx/N; %根据P109页(4.4.3),对rxx求期望生成自相关矩阵Rxx=rxx; %生成托布列兹阵(与本身共轭转置相等的方阵形式)V,S,U=svd(Rxx); %(U,S,V = svd (Rxx) %返回一个与X 同大小的对角矩阵S，两个酉矩阵U 和V，且满足Rxx= U*S*V。若Rxx为mn 阵，则U 为mm 阵，V为nn 阵。奇异值在S的对角线上，非负且按降序排列。这个分解得到信号自相关矩阵V,P109页4.4.4 S=V(:,1:2); %生成信号特征向量(主特征向量)P110页4.4.10,列从1到M G=V(:,3:N); %生成噪声特征向量,列从M到Nsyms w ; %定义符号 ,w(i)为2个正弦信号的频率ew=exp(-j.*(0:N-1).*(2*pi*w); %定义信号向量 见P112页上方e(w)定义Pw=1/(ew*(G*G)*ew); %求功率,见P112页公式(4.4.14)Pw=log10(abs(Pw); w=0:0.01:3; %频域范围figure;plot(w,subs(Pw); %频域图,subs为元素替换,没有矩阵替换后执行矩阵运算title(特征分解法music谱估计);grid on; 总结与展望 本文举