现代信号处理经典的功率谱估计_精品文档Word文档下载推荐.docx
《现代信号处理经典的功率谱估计_精品文档Word文档下载推荐.docx》由会员分享,可在线阅读,更多相关《现代信号处理经典的功率谱估计_精品文档Word文档下载推荐.docx(11页珍藏版)》请在冰豆网上搜索。
二、总体概述
本次实验分别使用经典的功率谱估计(如周期图法)与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;
%关闭之前的所有的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做向量的共轭乘积
figure
(2);
plot(w/pi/2,pw2);
%输出功率谱函数pw2图像
正弦信号x(n)自相关法功率谱估计'
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值下功率谱密度波形的分辨率高低。
2.2源代码
%AR模型的一般最小二乘估计
%-----------------------------------------
%清除workspace之前的变量
%关闭之前的图像
%清除命令行之前的文字
n=[1:
128];
%取定采样点n=1至128
%取定f1频率的值
%取定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
fori=1:
p
forj=1:
R(i,j)=m(q+i+j-1-p);
%构造一个p*p阶的自相关矩阵(Hankel矩阵)
%(课本P883.4.33a)
end
end
Rlegnth=size(R)%输出验证R矩阵的行列数的值
p%i=1~p
r(i)=m(q+i);
%定义一个1*p的向量,对应课本P883.4.22c
r=-r'
;
%对应课本P883.4.23Ra=-r
a=(inv(R'
*R)*R'
)*r;
%用LS方法求解a
a1=fliplr(a)%对应课本P883.4.22b,将a进行元素对调,使a1=[ap,...,a1]'
freqz(1,a1,128,1);
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求总体最小二乘解的算法步骤,编写程序。
(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算法)
%清除workspace之前的变量
%关闭之前的图像
%清除命令行之前的文字
%取定采样点n=1至128