傅里叶变换进行谱矩估计

### 使用傅里叶变换进行谱矩估计 #### 方法概述 频域分析中的一个重要工具是通过快速傅立叶变换(FFT)来计算信号的功率谱密度(Power Spectral Density, PSD)[^1]。对于给定的时间序列数据$x(t)$，其离散形式可以表示为$x[n]$。利用FFT算法可获得该时间序列对应的频率成分$X(f_k)=\sum_{n=0}^{N-1}{x[n]\cdot e^{-j2πf_kt_n}}$，其中$f_k=k/T_s(k=0,\cdots,N/2;T_s=NΔt)$代表各次谐波分量所对应的实际物理意义下的频率值[^2]。 基于上述得到的复数型频谱系数$|X(f)|^2/N$即反映了原始信号的能量分布情况——也就是常说的PSD函数图像；而所谓的“谱矩”，则是指对这一能量分布按照特定阶数求取统计特征的过程。具体来说： - 零阶矩（总功率）：$\int_0^\infty S_x(f)\mathrm{d}f=\frac{\Delta f}{N}\sum_{k=0}^{N/2}|X[k]|^2$ - 一阶矩（质心频率）：$\bar{f}= \left(\int_0^\infty fS_x(f)\mathrm{d}f\right)/\left(\int_0^\infty S_x(f)\mathrm{d}f\right)=(\Delta f/\sum |X[k]|^2 )\times (\sum k|X[k]|^2)$ - 二阶中心距（带宽）：$\sigma_f^2=(\Delta f)^2\times[\sum (k-\mu_f)^2|X[k]|^2]/\sum |X[k]|^2$ 以上各式中，$\Delta f=f_s/N$ 表示相邻两个样本点之间的间隔宽度，$f_s$ 是采样率，$N$ 则是指整个周期内采集到的数据总量[^3]。 #### Python 实现代码 下面给出一段Python代码用于展示如何应用numpy库完成从读入实测振动加速度记录直至最终得出所需各类谱矩指标全过程的操作方法: ```python import numpy as np from scipy.fft import fft def compute_spectrum_moments(signal, fs): """ 计算并返回零阶、一阶以及二阶谱矩 参数: signal : array_like 输入的一维时序信号. fs : float 采样频率(Hz). 返回: tuple of floats 包含三个元素分别对应于零阶矩(M0), 一阶矩(M1)， 和二阶矩(M2) """ N = len(signal) # 执行快速傅里叶变换获取双边幅度谱 X = fft(signal) # 取单边幅值平方作为功率谱密度估算 psd = ((np.abs(X[:N//2])**2)/(fs*N)) df = fs / N freqs = np.arange(N // 2)*df M0 = sum(psd * df) M1 = sum(freqs*psd * df) / M0 M2 = sum(((freqs-M1)**2)*psd * df) return M0, M1, M2 if __name__ == "__main__": t = np.linspace(0., 1., num=int(8e3)) sig = np.sin(2*np.pi*t)+0.5*np.random.randn(len(t)) # 构造测试用正弦波加上随机噪声混合而成的新号源 m0,m1,m2 = compute_spectrum_moments(sig, int(8e3)) print('Zeroth moment:',m0,'Hz') print('First moment or mean frequency:',m1,'Hz') print('Second central moment or bandwidth squared:',m2,'Hz²') ```