现代信号处理经典的功率谱估计
现代信号处理功率谱估计
现代信号处理功率谱估计
式中, p(x)是X的概率密度函数，对于离散随机序列, 概率密度函 数用联合概率密度函数代替。显然,熵代表一种不确定性， 最大 熵代表最大的不确定性， 或者说最大的随机性。下面我们研究 对于有限的自相关函数值不作任何改变，对于未知自相关函数 用最大熵原则外推，即不作任何附加条件的外推方法。 假设x(n) 是零均值正态分布的平稳随机序列，它的N维高斯概率密度函数 为 p ( x 1 ,x 2 , ,x N ) ( 2 π ) N /2 (d R x( N x e )1 /2 ) e t x 1 2 X H p ( R x( N x) 1 X )
rxx(1)
rxx(2)
rxx(0) rxx(1)
rxx(N
1)
rxx(N
2)
0
rxx(N1) rxx(N)
rxx(1)
可以看出AR模型得到的结果与按最大熵外推rxx(N+1)得到的结果 一致，这就证明了当x(n)为高斯分布时的最大熵谱估计与AR模型
法是等价的。
上式（4.6.8)是rxx(N+1)的一次函数，由此可解得rxx(N+1)。再 用类似的方法求得rxx(N+2)， rxx(N+3)，┄，然后确定功率谱估计。
式中det(Rxx(N))表示矩阵Rxx(N)的行列式，由上式表明为使熵最 大，要求det(Rxx(N)最大。
现代信号处理功率谱估计
若已知N+1个自相关函数值rxx(0),rxx(1),…,rxx(N)，下面用最 大熵方法外推rxx(N+1)。设rxx(N+1)确实是信号自相关函数的第 N+2个值，根据自相关函数的性质，由N+2个自相关函数组成 的矩阵为
经典功率谱估计
雷达和声呐系统
目标检测
在雷达和声呐系统中，经典功率谱估计常被用于目标检测。通过对接收到的信号进行功率 谱分析，可以判断是否存在目标以及目标的位置和速度等信息。
距离和速度测量
在雷达和声呐系统中，经典功率谱估计还可以用于距离和速度测量。通过对接收到的信号 进行功率谱分析，可以估计出目标与系统之间的距离和相对速度。
信号分类
在雷达和声呐系统中，经典功率谱估计还可以用于信号分类。通过对接收到的信号进行功 率谱分析，可以判断目标的类型，例如区分飞机、船舶或车辆等不同类型目标。
05 经典功率谱估计的改进方 法
基于小波变换的功率谱估计
1
小波变换能够将信号分解成不同频率和时间尺度 的分量，从而更好地揭示信号的内在结构和特征。
然而，这些方法通常需要较长 的数据长度和较为复杂的计算 过程，对于短数据和实时处理 的应用场景具有一定的局限性 。
研究展望
01
随着信号处理技术的发展，经典功率谱估计方法仍有进一步优化的空 间。
02
针对短数据和实时处理的应用场景，研究更为快速、准确的功率谱估 计方法具有重要的实际意义。
03
结合机器学习和人工智能技术，探索基于数据驱动的功率谱估计方法 是一个值得关注的方向。
优点
能够提供较高的频率分辨率和较低的估计误差。
原理
格莱姆-梅尔谱估计利用了信号的模型参数，通过 构造一个模型函数来描述信号的频率响应特性， 并求解该函数的极值问题得到信号的功率谱。
缺点
需要预先设定模型函数的形式和参数，且计算复 杂度较高。
03 经典功率谱估计的优缺点
优点
01
02
03
算法成熟
经典功率谱估计方法经过 多年的研究和发展，已经 相当成熟，具有较高的稳 定性和可靠性。
现代信号处理 总结1
第1章 离散时间信号与系统1、 傅里叶分析和Z 变换的区别、缺陷、特点关系：点数为N 的有限长序列x(n)的Z 变换为X(z),而其离散傅里叶变换为X(k)，两者均表示了同一有限长序列x(n)的变换,它们之间的关系是:对z 变换在单位圆上取样可得DFT 。
而DFT 的内插就是变换。
傅里叶变换优缺点(1) 傅里叶变换缺乏时间和频率的定位功能 (2) 傅里叶变换对于非平稳信号的局限性(3) 傅里叶变换在时间和频率分辨率上的局限性傅立叶变换是最基本得变换，由傅里叶级数推导出。
傅立叶级数只适用于周期信号，把非周期信号看成周期T 趋于无穷的周期信号，就推导出傅里叶变换，能很好的处理非周期信号的频谱。
但是傅立叶变换的弱点是必须原信号必须绝对可积，因此适用范围不广。
Z 变换的本质是离散时间傅里叶变换（DTFT ），如果说拉普拉斯变换专门分析模拟信号，那Z 变换就是专门分析数字信号，Z 变换可以把离散卷积变成多项式乘法，对离散数字系统能发挥很好的作用。
Z 变换看系统频率响应，就是令Z 在复频域的单位圆上跑一圈，即Z=e^(j2πf)，即可得到频率响应。
2、系统的记忆性、因果性、可逆性(1)记忆性如果系统在任意时刻n0的响应仅与该时刻的输入f(n0)有关，而与其它时刻的输入无关，则称该系统为非记忆系统（或系统无记忆性），否则称为记忆系统。
系统的记忆性有时也被称为动态特性。
该特性强调系统的响应是否仅与当前时刻的输入有关。
对于无记忆LTI 系统，其系统冲激响应为，其中()()h n K n δ=，K 为一常数。
由于系统频率响应是冲激响应的傅氏变换、系统函数为系统冲激响应的z 变换，因此，无记忆LTI 系统的系统频率响应和系统函数分别为H(ω)=K ，H(z)=K 。
(2) 因果性如果系统任意时刻的响应与以后的输入无关，则该系统称为因果系统（或系统具有因果性），否则为非因果系统。
该特性强调的是，系统的响应是否与未来的输入有关。
功率谱估计的方法
功率谱估计的方法
功率谱估计是信号处理中常用的一种方法，用于分析信号在频域内的特点，通常可以分为以下几种方法：
一、经典方法
1.傅里叶变换法：将时域信号通过傅里叶变换变换到频域，然后计算功率谱密度。
2.自相关法：通过自相关函数反映信号的统计平稳性，然后通过傅里叶变换计算功率谱密度。
3.周期图法：将信号分解为若干个周期波形，然后对每个周期波形进行傅里叶变换计算周期功率谱，最后汇总得到整个信号的功率谱。
二、非经典方法
1. 时-频分析法：如短时傅里叶变换（STFT）、小波变换等，将信号分解为时域和频域两个维度的分量，从而可以分析信号在时间和频率上的变化。
2. 基于协方差矩阵的特征值分解法：通过建立协方差矩阵，在张成空
间中求解特征向量，从而达到计算信号功率谱的目的。
3. 基于频率锁定法：如MUSIC法、ESPRIT法等，是一种利用特定信号空间中的特定模式进行处理的方法。
以上方法各有特点，根据实际需求选择不同的方法可以得到相应的功率谱估计结果。
(周围)现代信号处理基础04-现代功率谱估计(上)
• ARMA模型 H (z) B(z) A( z )
参数法谱估计的理论基础
重庆邮电大学通信学院
谱分解定理的推论
任何平稳随机信号x(n)都可以看成由白噪声序列 {u(n)} 激励 一个因果和稳定的线性时不变系统H(z)产生的输出。
任何有限方差的平稳ARMA过程可以分为完全随机的部分 和确定的部分，对应的功率谱为连续的和离散的冲激信号。
FT
FT () fT (t)
IFT
fT
(t )
f
(t), T 2
t
T 2
0,
其他
• 平稳离散随机信号x(n)的自相关函数与功率谱密度之间为一 对傅立叶变换
物理意义：功率Rxx (0)
功率在ω上的分布
重庆邮电大学通信学院
如果随机信号是各态遍历的，相关函数可以由一个取样时间 序列用时间平均来取代统计平均。
AR模型法功率谱估计：
AR(p)模型的Yule - Walker方程组：
R (0) R (1)
R (1)
R (0)
R (2) R (1)
R (p) R (p-1)
R (2) R (1) R (0)
R (p-2)
R (p) 1 2
R
(p-1)
重庆邮电大学通信学院
AR谱估计的性质2：与最大熵谱估计等效
结论：
AR谱估计相当于对自相关函数以最大熵为原则进行 外推后进行傅立叶变换的结果。
AR谱估计相当于在p+1个自相关函数值确定的情况 下，以功率谱密度最平坦为准则得到的估计结果。
重庆邮电大学通信学院
生医信号处理课件6经典功率谱
（2）方差 Var[S (e j )] Var[S (e j )] S 2 (e j ) M per x 频率分辨率及旁瓣泄漏
Re s[SM (e j )] ()3dB
二、平均法（Bartlett法）
Bartlett提出，将 x(n)分为长L，互不重叠的k段子序列，N=kL,
3、讨论
m M 2 w ( m) M
M↑偏差↓
M↑方差↑ 一般M=N/5
五、谱估计技术的性能指标
1、变异性γ （归一化方差）
2、品质因数μ
4 2 x
j
j
结论：非一致估计
例
三、周期图的随机起伏
Cov[S per (e
j1
), S per (e
j2
sin[(k l ) ] 2 )] [ ] (k l ) N sin[ ] N
4 x
一、数据加窗（修正周期图）
数据窗
6.4.3 功率谱估计的改进
E[ S per (e )]
j
j
1 E[ S BT (e )] S x ( e j ) W ( e j ) 2
j
N ,
W (e j ) ()
是渐近无偏
2、方差
if 1 M N
j
1 2 Var[ S BT (e )] S x (e j ) W 2 (e j ( u ) )du 2N 1 j S x (e ) N
Bartlett法和Welch法分别对周期图和修正周期图进行平均， 从而达到减少方差的目的。 Blackman-Tukey法为了减少周
期图的方差，对自相关序列的估计进行加窗处理，从而减
少自相关序列的估计中那些不可靠的估计值对周期图的贡 献。
功率谱估计的经典方法PPT课件
吉林大学通信工程学院信息科学实验室
6
时间平均
（11）一个平稳随机过程的一个取样序列的时间平均等于它的集合平
均，则称它是遍历性随机过程。时间平均记为 x(n) ，则取样序列的算术
平均值和时间取样自相关序列定义为
x(n) lim 1
功率谱估计的经典方法
版权所有
吉林大学通信工程学院信息科学实验室
1
离散随机过程
为了描述随机变量，引入了概率分布函数、概率密度函数以及随机变 量的数字特征。这些函数或参数都是针对一维随机变量定义的。统称一 维统计特征。
但对于离散随机过程，因为它是由无限多个随机变量构成的时间序列
xn, n ，因此为完整地描述它，仅知道随机变量的特征是不
Syy(z)
Ryy(m) zm
Rxx(m
p)Rhh (
p)
zm
m
m p
Rxx(n)Rhh ( p)
z n z p
Sxx(z)Shh (z)
m n
S
xx
(
z
)H
(
z
)
H
(
z
1
)
协方差序列的z变换
Sxx(z) Cxx(m) zm , m
称为平稳随机过程的功率谱。在今后的讨论中总假设随机信号的均值为
零，所以有
Sxx(z) Rxx(m) zm , m
由于 Rxx(m) Rxx(m) ，则有 Sxx (z) Sxx (z 1) 。
经典功率谱估计和现代功率谱估计
三、实验原理
直 接 法 （ 周 期 图 法 ） ： 它 将 信 号 X(n) 的 N 点 样 本 XN ( ejw ) 为 XN n 的傅氏变换，取其幅值平方除以 N 作为 x(n)真实的功率谱 P( ejw ) 的估计。 P(ejw )表示。 间接法（自相关法、 BT 法）：先由 XN n 估计出自相关函数 r m 。对 r m 求傅氏 变换得到 XN n 的功率谱： PBT (ω)这是对 P( ω) 的估计。
相关图形：
相关函数法 35
30
25
20
15
10
5
0
50
100
150
200
250
300
周期图法 50
0
-50
-100
0
50
100
150
200
250
300
50
0
-50
-100
0
50100Fra bibliotek150
200
250
300
周期图法功率谱估计 20 10
Burg Power Spectral Density Estimate 60
40 0
功 率 谱 密 度 （ dB/Hz ） Power/frequency (dB/Hz)
20
-10 -20 -30 -40
0
-20
-40 -50 -60 -60
0
50 频 率 （ Hz ）
100
0
50 Frequency (Hz)
100
六、实验总结 可以看出直接法与间接法的方差性能都比较差， 为了追求谱线平滑， 就要以牺牲分辨率为代价的。
ylabel('¹¦ÂÊÆ×Ãܶȣ¨dB/Hz£©') title('ÖÜÆÚͼ·¨¹¦ÂÊÆ×¹À¼Æ') order1=50; range='half'; magunits='dB'; subplot(1,2,2) pburg(xn,order1,nfft,Fs,range);%burgËã·¨µÃ¹¦ÂÊÆ×
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
《现代信号处理》姓名：李建强学号：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; %关闭之前的所有的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.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至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>=pfor i=1:pfor j=1:pR(i,j)=m(q+i+j-1-p); %构造一个p*p阶的自相关矩阵（Hankel矩阵）%(课本P88 3.4.33a)endendRlegnth=size(R) %输出验证R矩阵的行列数的值for i=1:p % i=1~pr(i)=m(q+i); %定义一个1*p的向量，对应课本P88 3.4.22cendr=-r'; %对应课本P88 3.4.23 Ra=-ra=(inv(R'*R)*R')*r; %用LS方法求解aa1=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至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>=pfor i=1:pfor j=1:pR(i,j)=m(q+i+j-1-p); %构造一个pxp阶的自相关矩阵（Hankel矩阵）%(课本P88 3.4.33a)endendRlegnth=size(R) %输出验证R矩阵的行列数的值for i=1:p %i=1~pr(i)=m(q+i); %定义一个1*p的向量，对应课本P88 3.4.22cendB=[-r',R]; %对应P94 3.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，对应课本P95 3.4.55for j=1:pfor i=1:p+1-Pdjj=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，构造vijS= S+djj*vij*vij'; %对应课本P96 3.4.56，计算矩阵S的二重级数求和endendSni=inv(S); %对应课本P96 算法3.4.1步骤4，求S逆矩阵a=zeros(1,P); %对应课本P88 3.4.22b，构造a矩阵for i=1:Pa(1,i)=Sni(i+1,1)/Sni(1,1); %对应课本P96 3.4.59，求出矩阵a=[a1,...,ap]' enda1=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、、周期图法的分辨率低，不能适应高分辨率功率谱估计的需要，与之相比，参数化谱估计可以提供比周期图高得多的频率分辨率。