现代信号处理经典的功率谱估计

。-可编辑修改-《现代信号处理》姓名：李建强学号：201512172087专业：电子科学与技术作业内容：在MATLAB平台上对一个特定的平稳随机信号进行经典功率谱估计和现代功率谱估计的比较一、前言功率谱估计是信息学科中的研究热点,在过去的30多年里取得了飞速的发展。在许多工程应用中，它能给出被分析对象的能量随频率的分布情况。平滑周期图是一种计算简单的经典方法，它的主要特点是与任何模型参数无关，但估计出来的功率谱很难与信号的真是功率谱相匹配。与周期图方法不同，现代谱估计主要是针对经典谱估计(周期图和自相关法)的分辨率低和方差性能不好的问题而提出的。其使用参数化的模型，能够给出比周期图方法高得多的频率分辨率。其内容极其丰富，涉及的学科和领域也相当广泛，按是否有参数大致可分为参数模型估计和非参数模型估计，前者有AR模型、MA模型、ARMA模型、PRONY指数模型等;后者有最小方差方法、多分量的MUSIC方法等。二、总体概述本次实验分别使用经典的功率谱估计（如周期图法）与AR模型法对某一特定的平稳随机信号进行其功率谱估计，由图像得到信号的频率。利用MATLAB平台，直观形象地观察并比较二者估计效果的区别，以便于加深对功率谱估计的理解和掌握。三、具体的实现步骤1、经典法功率谱估计周期图法又称直接法，它是从随机信号x(n)中截取N长的一段，把它视为能量有限的。-可编辑修改-真实功率谱的估计的一个抽样。1.1、实现步骤（1）、模拟系统输出参数x(n)=A*sin(2πf1*n)+B*sin(2πf2*n)，包括序列长度N（128或512或1024，加性高斯白噪声（AGWN）功率一定，设置A，B，f1，f2，n的值。（2）、应用周期图法（不加窗）对信号的功率谱密度进行估计，使用直接法在MATLAB平台上进行编程实现。（3）、输出相应波形图，进行观察，记录。1.2MATLAB源代码实现clearall;%清除工作空间所有之前的变量closeall;%关闭之前的所有的figureclc;%清除命令行之前所有的文字n=1:1:128;%设定采样点n=1-128f1=0.2;%设定f1频率的值0.2f2=0.213;%设定f2频率的值0.213A=1;%取定第一个正弦函数的振幅B=1;%取定第一个正弦函数的振幅a=0;%设定相位为0x1=A*sin(2*pi*f1*n+a)+B*sin(2*pi*f2*n+a);%定义x1函数，不添加高斯白噪声x2=awgn(x1,3);%在x1基础上添加加性高斯白噪声，信噪比为3,定义x2函数temp=0;%定义临时值，并规定初始值为0temp=fft(x2,128);%对x2做快速傅里叶变换pw1=abs(temp).^2/128;%对temp做经典功率估计。-可编辑修改-k=0:length(temp)-1;w=2*pi*k/128;figure(1);%输出x1函数图像plot(w/pi/2,pw1)%输出功率谱函数pw1图像xlabel('信号频率/Hz');ylabel('PSD/傅立叶功率谱估计');title('正弦信号x(n)添加高斯白噪声后的,周期图法功率频谱分析');grid;%-------------------------------------------------------------------------pw2=temp.*conj(temp)/128;%对temp做向量的共轭乘积k=0:length(temp)-1;w=2*pi*k/128;figure(2);plot(w/pi/2,pw2);%输出功率谱函数pw2图像xlabel('信号频率/Hz');ylabel('PSD/傅立叶功率谱估计');title('正弦信号x（n）自相关法功率谱估计');grid;1.3matlab仿真图形（1）、用直接法，功率谱图像,采样点N=128。。-可编辑修改-（2）用直接法，功率谱图像,采样点N=512。1.4、经典功率谱估计分析。-可编辑修改-当采样的点数为N=128时，此时采样的得到的图像分辨力很低，并且分辨率也比较低，这就导致了功率谱图像只能看到一个峰值点。采样点数为N=512时，此时，分辨力和分辨率比较高，可以清楚的区分到两个峰值点的横坐标，此时的横坐标就是信号的频率。但是这是以牺牲效率为代价的，采样的点数越多，所花的时间越长，这在实际的工程中是不切合实际的，因此，在我们估计随机信号的频率的时候，要合理的采取样本点数，尽可能的采取多的样点，来接近真实的信号频率，也要考虑实际的效率问题。2、AR模型一般最小二乘法谱分析方法要求ARMA模型的阶数和参数以及噪声的方差已知.然而这类要求在实际中是不可能提供的,即除了一组样本值x(1),x(2),…,x(T)以供利用(有时会有一定的先验知识)外,再没有其它可用的数据.因此必须估计有关的阶数和参数,以便获得谱密度的估计。2.1实现步骤（1）、模拟系统输出参数y=A*sin(2πf1*n)+B*sin(2πf2*n)，包括序列长度N,加性高斯白噪声（AGWN）,设置A，B，f1，f2，n的值。（2）、应用AR模型一般最小二乘法对信号进行功率谱估计，编写程序。取定|B(z)|=1，构造AR模型，然后不断变换p的值，观察不同p值下功率谱密度波形的分辨率高低。（3）、输出相应波形图，进行观察，记录。2.2源代码%AR模型的一般最小二乘估计%-----------------------------------------clearall;%清除workspace之前的变量closeall;%关闭之前的图像。-可编辑修改-clc;%清除命令行之前的文字n=[1:128];%取定采样点n=1至128f1=0.2;%取定f1频率的值f2=0.213;%取定f2频率的值(根据f1与f2之差=2*pi/n=0.0491）A=sqrt(20);%取定第一个正弦函数的振幅B=sqrt(2);%取定第一个正弦函数的振幅x=A*sin(2*pi*f1*n)+B*sin(2*pi*f2*n);%定义x函数noise=0+1*randn(1,length(n));%添加均值为0、方差为1的高斯白噪声xn=x+noise;%在x1基础上添加加性高斯白噪声，定义xn函数m=xcorr(xn);%m为xn的自相关函数（序列)%-----------------------------------------p=100;%取定R的阶数,更改p的值，观察相对应的谱估计q=125;%此处一定要满足q=pfori=1:pforj=1:pR(i,j)=m(q+i+j-1-p);%构造一个p*p阶的自相关矩阵（Hankel矩阵）%(课本P883.4.33a)end。-可编辑修改-endRlegnth=size(R)%输出验证R矩阵的行列数的值fori=1:p%i=1~pr(i)=m(q+i);%定义一个1*p的向量，对应课本P883.4.22cendr=-r';%对应课本P883.4.23Ra=-ra=(inv(R'*R)*R')*r;%用LS方法求解aa1=fliplr(a)%对应课本P883.4.22b，将a进行元素对调,使a1=[ap,...,a1]'figure(1);freqz(1,a1,128,1);title('AR模型的一般最小二乘估计');legend(strcat('AR阶数=',int2str(p)));gridon;2.3matlab仿真图形。-可编辑修改-（1）、当p=4时，信号的功率谱密度波形：（2）、当p=64时，信号的功率谱密度波形：。-可编辑修改-（3）、当p=100时，信号的功率谱密度波形：3、AR模型的总体最小二乘法。-可编辑修改-一般最小二乘法会带来两个问题：其一，必须重新列出方程组，使它只包含p个未知数；其二，求解Ax=b的最小二乘方法只认为b含有误差，但实际上系数矩阵A也含有误差。因此，引入总体最小二乘法（SVD-TLS）可以比一般最小二乘法更合理地同时考虑A和b的误差或扰动。3.1实现步骤（1）、模拟系统输出参数y=Asin(2πf1*n)+Bsin(2πf2*n)，包括序列长度N,加性高斯白噪声（AGWN）,设置A,B,f1,f2,n的值。（2）、应用AR模型的总体最小二乘法对信号进行功率谱估计，按照课本P96求总体最小二乘解的算法步骤，编写程序。取定|B(z)|=1，构造AR模型，然后不断变换p的值，观察不同p值下功率谱密度波形的分辨率高低。（3）、输出相应波形图，进行观察，记录。（4）、算法步骤（SVD-TLS算法）步骤1：计算增广矩阵B的SVD，并存储奇异值和矩阵V；步骤2：确定增广矩阵B的有效秩p；步骤3：利用式（3.4.56）和式（3.4.53）计算矩阵S（p）；步骤4：求S（p）的逆矩阵S-(p)，并由式（3.4.59）计算未知参数的总体最小二乘估计。3.2源代码%实现AR模型的总体最小二乘估计（SVD-TLS算法）%-----------------------------------------clearall;%清除workspace之前的变量closeall;%关闭之前的图像。-可编辑修改-clc;%清除命令行之前的文字n=[1:128];%取定采样点n=1至128f1=0.2;%取定f1频率的值f2=0.213;%取定f2频率的值A=sqrt(20);%取定第一个正弦函数的振幅B=sqrt(2);%取定第二个正弦函数的振幅x=A*sin(2*pi*f1*n)+B*sin(2*pi*f2*n);%定义x函数noise=0+1*randn(1,length(n));%添加均值为0、方差为1的高斯白噪声xn=x+noise;%在x1基础上添加加性高斯白噪声，定义xn函数m=xcorr(xn);%m为xn的自相关函数（序列)%-----------------------------------------p=100;%取定R的阶数,更改p=4，64，100，的值，观察%相对应的谱估计q=125;%此处一定要满足q=pfori=1:pforj=1:pR(i,j)=m(q+i+j-1-p);%构造一个pxp阶的自相关矩阵（Hankel矩阵）%(课本P883.4.33a)end。-可编辑修改-endRlegnth=size(R)%输出验证R矩阵的行列数的值fori=1:p%i=1~pr(i)=m(q+i);%定义一个1*p的向量，对应课本P883.4.22cendB=[-r',R];%对应P943.4.45b中的B[U,K,V]=svd(B);%由P96算法3.4.1步骤1求得增广矩阵B的%SVD，并存储奇异值和矩阵VP=rank(B);%由P96算法3.4.1步骤2求得增广矩阵B的有%效秩，定义为PS=zeros(P+1);%构造一个(p+1)*(p+1)维的矩阵S，对应课本P953.4.55forj=1:pfori=1:p+1-Pdjj=K(j,j)*K(j,j);%对应课本P963.4.56，构造djj,并求其平方vij=V(i:i+p,j);%对应课本P963.4.56和课本P953.4.53，构造vijS=S+djj*vij*vij';%对应课本P963.4.56，计算矩阵S的二重级数求和endend。-可编辑修改-Sni=inv(S);%对应课本P96算法3.4.1步骤4，求S逆矩阵a=zeros(1,P);%对应课本P883.4.22b，构造a矩阵fori=1:Pa(1,i)=Sni(i+1,1)/Sni(1,1);%对应课本P963.4.59，求出矩阵a=[a1,...,ap]'enda1=fliplr(a)%对应课本P883.4.22b，将a进行元素对调,使a1=[ap,...,a1]'figure(1);freqz(1,a1,128,1);%求出信号的幅频响应和相频响应波形title('AR模型的总体最小二乘（SVD-TLS）估计');legend(strcat('AR（p）阶数p='