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

《现代信号处理》

姓名：李建强

学号：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.2 MATLAB源代码实现

clear all; %清除工作空间所有之前的变量

close all; %关闭之前的所有的figure

clc; %清除命令行之前所有的文字

n=1:1:128; %设定采样点n=1-128

f1=0.2; %设定f1频率的值0.2

f2=0.213; %设定f2频率的值0.213

A=1; %取定第一个正弦函数的振幅

B=1; %取定第一个正弦函数的振幅

a=0; %设定相位为0

x1=A*sin(2*pi*f1*n+a)+B*sin(2*pi*f2*n+a); %定义x1函数，不添加高斯白噪声x2=awgn(x1,3); %在x1基础上添加加性高斯白噪声，信噪比为3,定义x2函数temp=0; %定义临时值，并规定初始值为0

temp=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.3 matlab仿真图形

（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模型的一般最小二乘估计

%-----------------------------------------

clear all; %清除workspace之前的变量

close all; %关闭之前的图像

clc; %清除命令行之前的文字

n=[1:128]; %取定采样点n=1至128

f1=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>=p

for i=1:p

for j=1:p

R(i,j)=m(q+i+j-1-p); %构造一个p*p阶的自相关矩阵（Hankel矩

阵）

%(课本P88 3.4.33a)

end

end

Rlegnth=size(R) %输出验证R矩阵的行列数的值

for i=1:p % i=1~p

r(i)=m(q+i); %定义一个1*p的向量，对应课本P88 3.4.22c

end

r=-r'; %对应课本P88 3.4.23 Ra=-r

a=(inv(R'*R)*R')*r; %用LS方法求解a

a1=fliplr(a) %对应课本P88 3.4.22b，将a进行元素对调,使a1=[ap,...,a1]'

figure(1);

freqz(1,a1,128,1);

title('AR模型的一般最小二乘估计');

legend(strcat('AR阶数=',int2str(p)));

grid on;

2.3 matlab仿真图形

（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算法）

%-----------------------------------------

clear all; %清除workspace之前的变量

close all; %关闭之前的图像

clc; %清除命令行之前的文字

n=[1:128]; %取定采样点n=1至128

f1=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>=p

for i=1:p

for j=1:p

R(i,j)=m(q+i+j-1-p); %构造一个pxp阶的自相关矩阵（Hankel矩阵）

%(课本P88 3.4.33a)

end

end

Rlegnth=size(R) %输出验证R矩阵的行列数的值

for i=1:p %i=1~p

r(i)=m(q+i); %定义一个1*p的向量，对应课本P88 3.4.22c

end

B=[-r',R]; %对应P94 3.4.45b中的B

[U,K,V]=svd(B); %由P96 算法3.4.1步骤1 求得增广矩阵B的%SVD，并存

储奇异值和矩阵V

P=rank(B); %由P96 算法3.4.1步骤2 求得增广矩阵B的有%效秩，定

义为P

S=zeros(P+1); %构造一个(p+1)*(p+1)维的矩阵S，对应课本P95 3.4.55

for j=1:p

for i=1:p+1-P

djj=K(j,j)*K(j,j); %对应课本P96 3.4.56，构造djj,并求其平方

vij=V(i:i+p,j); %对应课本P96 3.4.56和课本P95 3.4.53，构造vij

S= S+djj*vij*vij'; %对应课本P96 3.4.56，计算矩阵S的二重级数求和

end

end

Sni=inv(S); %对应课本P96 算法 3.4.1步骤4，求S逆矩阵

a=zeros(1,P); %对应课本P88 3.4.22b，构造a矩阵

for i=1:P

a(1,i)=Sni(i+1,1)/Sni(1,1); %对应课本P96 3.4.59，求出矩阵a=[a1,...,ap]'

end

a1=fliplr(a) %对应课本P88 3.4.22b，将a进行元素对调,使a1=[ap,...,a1]'

figure(1);

freqz(1,a1,128,1); %求出信号的幅频响应和相频响应波形

title('AR模型的总体最小二乘（SVD-TLS）估计');

legend(strcat('AR（p）阶数p=',int2str(p)));

grid on;

3.3、matlab仿真图形

（1）、当p=4时，信号的功率谱波形

（2）、当p=64时，信号的功率谱波形

（3）、当p=100时，信号的功率谱波形

四、总结：

以上分别用了周期图法、AR模型的参数化估计（LS算法和SVD-TLS）算法对同一观测信号进行了功率谱的估计，通过仿真结果的对比，可以得出以下结论：

1、、周期图法的分辨率低，不能适应高分辨率功率谱估计的需要，与之相比，参数化谱估计可以提供比周期图高得多的频率分辨率。应用周期图法不能够分辨出f1和f2两处功率谱的尖峰，应用最小二乘法则可以明显的分辨出f1=0.2和f2=0.213两处的功率谱，分辨率高。

2、LS算法和SVD-TLS算法比较，仿真波形的误差较大，这是由于LS算法的ARMA的阶数P是任意设定的，并没有遵循严格的理论依据，而且在Ax=b中，LS算法只是考虑了b的误差和扰动，而SVD-TLS算法则是综合考虑了A和b的误差与扰动。所以LS算法存在较大误差，而SVD-TLS算法则有较好的结果。

3、无论是LS算法还是SVD-TLS算法，阶数R(M,p,q)的设置不同，会导致得出不同效果的仿真波形。比如当p=4时（固定q=125），波形的尖峰并没有明显地区分开来，但当p=64和p=100时，就会看到尖峰分离地非常明显了。所以，在LS和 SVD-TLS算法中，M，p,q（q>=p）的设置应尽量大一些，这样才能得到分辨率良好的功率谱波形。