经典法谱估计（1）周期图法MATLAB及Python实现

原理

对于一个无限长序列x(n)，其傅里叶变换DTFT（也可称作离散时间傅里叶变换）及其反变换的定义式为：

但是在我们的分析处理过程中，只能分析处理离散化的信号，也即N点取样，对于一个长度为N的有限长序列x(n)，其离散傅里叶变换及其反变换的定义式为：

或者也可以写为

式中x(n)和X(k)是一个有限长序列的离散傅里叶变换对。长度为N的有限长序列x(n)，其离散傅里叶变换X(k)是一个有限长频域序列，其长度依然为N。对X(k)的表达式分析，可知其周期为N，我们现仅是截取一个周期内的频域信号进行分析。由于序列x(n)的离散傅里叶变换X(k)具有周期性，因而这种功率谱也具有周期性，所以称为周期图法。
所谓功率谱密度，即为在单位谱上的功率分布，因此，我们可以推知

即得到原序列x(n)的功率谱。
这种方法以很简洁的方式求得了功率谱，但是其缺点在于：只能用于求解确定性信号的功率谱密度，而不能用于求解随机信号的功率谱密度。原因在于：随机序列x(n)在Z域不满足绝对可和条件，也就不存在其DFT。例如给一个确定性信号加入了高斯白噪声，原信号序列是绝对可和的，故能够用周期图法，但是加入噪声后的新的信号序列，它的DFT不存在，如果只是单纯不考虑这个而取N个值来求傅里叶变换，所得到的功率谱每次是不一样的，已经失去了其意义，因此周期图法不再适用，而要选择后面提到的相关函数等方法。

程序和结果

（出于维护版权原因，此处只放程序截图）
MATLAB：
程序：


结果：

Python
程序：

结果：

分析

由上图可见，我给程序输入的N为512，此时我们可以从功率谱图中看到四个峰，呈现对称分布，信号中f1=0.1，f2=0.13，在图中我们可以看到对应的位置出现了峰值，由于给定幅度大小不同，故峰值大小有很大差别。且由于输入的x(n)为一个实序列，根据“实偶虚奇”，它所对应的DFT为偶对称，这也就是为什么会出现左右对称的四个峰。由于只有确定性信号有DFT，随机信号不满足绝对可和而导致其DFT不存在，故我在这里没有给信号加入噪声，而是求的原信号的功率谱，没有考虑噪声影响，峰值位置理论上是在两个频率处。

（因原博客是以word形式写的，CSDN不支持Mathtype公式，故部分地方直接采用了截图形式）
下载链接：
经典法功率谱估计及MSE随SNR曲线