经典法谱估计（3）窗函数法Python实现

原理

我们在做FFT变换过程中只能对有限长度的时域数据进行变换，因此，需要对时域信号进行信号截断。即使是周期信号，如果截断的时间长度不是周期的整数倍（周期截断），那么，截取后的信号将会存在泄漏。为了将这个泄漏误差减少到最小程度，我们需要使用加权函数，也叫窗函数。加窗主要是为了使时域信号似乎更好地满足FFT处理的周期性要求，减少泄漏。加窗实质是用一个所谓的窗函数与原始的时域信号作乘积的过程（频域卷积）使得相乘后的信号似乎更好地满足傅立叶变换的周期性要求。例如，一个窗函数可以定义为：

g(n)是窗函数，N是窗函数的长度。待分析的数据x(t)则表示为：

如此处理之后，再利用上述方法二，求x(n)的自相关函数再求自相关函数的傅里叶变换从来得到功率谱，相对于原先的结果来看，泄露会有所减少。
几种常见的窗函数：


我们可以调用MATLAB中的wvtool函数用来显示窗的形状和频域图形，以及得到其对应的参数。
例如：哈明窗：wvtool(hamming(64));

汉宁窗：wvtool(hann(64));

程序及结果

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



结果：

分析

由上图可见，我给程序输入的N为512。在上一种方法的基础上，对序列进行了加窗处理，我使用了汉宁窗，由功率谱图比较可以看出加窗后主瓣变窄，但同时它使得幅值失真，需要后期对频谱进行修正。修正分幅值修正和能量修正，如果是单条谱线则为幅值修正；如果是宽带则为能量修正。然后我继续使用上一种方法中提到的寻峰方式，在Python中没有已经封装好的库函数来寻峰，故我使用了分段求最大值的方式。由图中可以知道，峰值在0.1和0.13附近，那么就可以分别寻找[0,0.11],[0.11,0.5]范围内的峰值对应的索引值再转化为频率即可。在噪声功率远小于信号功率的情况下，在很接近信号频率的附近，噪声产生的扰动应该较小，几乎可以忽略，所以在信号频率附近的峰值应该可以说几乎是纯信号的峰值，而没有噪声引起的小峰，可想而知，区间内的峰值应该是信号的频率所在位置。由图可知，加窗前后，对于频率的估计值并没有太大的影响，这是因为取的N值较大，本身不加窗的时候也能有较好的功率谱估计性能。如果把N值取的小一点，例如取N=128时，有：

比较可以看出加窗与不加窗的区别，可以看出加窗的优越性。

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