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

《现代信号处理》

姓名：李建强

学号：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矩阵的行列数的值