傅里叶变换和功率谱密度

一、离散傅里叶变换

假设有信号$x(t)$，我们需要分析它的频率成分。

首先对$x(t)$采样，得到N个离散的$x_n$。做离散傅里叶变换，得到频域信号$X(f)$的采样

$$X_k=\frac{T}{N}\sum^{N-1}_{n=0}x_n \exp(-2\pi ikn/N).$$

(下标和自变量的关系是$t=n/N T，f=k/T, T$ 为采样时长。)

由上式可以看出，$X_k=X_{k+N}=-X_{k+N/2}$，因此可以减少计算量，使用快速傅里叶方法。此外，$x(t)$通常是实数，所以$X_k=X_{N-k}^{*}$，得到的频谱图有一定对称性，我们只取一半的信息，绘制单边谱：

import numpy as np
from scipy.fftpack import fft
import matplotlib.pyplot as plt
np.random.seed(2)

N = 128 # 采样点数
T=1 # 信号长度；信号长度整除周期时不会发生信号泄漏
dt=T/N # 采样间隔
Fs=N/T # 采样频率, 须大于两倍的最高频率 
t=np.arange(0,T,dt)

signal=5+2*np.sin(2*np.pi*20*t)+3*np.sin(2*np.pi*30*t)+4*np.sin(2*np.pi*40*t)+np.random.randn(len(t))  # 采集的信号
fft_data = fft(signal)*dt # 快速傅立叶变换
fft_amp0 = np.array(np.abs(fft_data)*2) 
fft_amp0[0]=0.5*fft_amp0[0] # 0频率的幅值除以点数N，非零频率除以 N/2
list0 = np.array(range(0, N))
freq0 = list0/T # 频率坐标
plt.plot(freq0, fft_amp0)
plt.xlim(0,int(Fs/2))
plt.xlabel('frequency (Hz)')
plt.ylabel(' Amplitude (unit/Hz)')
plt.show()

由上图可以看出，$x(t)$由5个信号叠加，频率为0，20Hz，30Hz和40Hz, 将对应纵轴的值除以信号长度$T$，得到这五个信号在时域中的幅值为5，2，3，4。

二、功率谱密度

功率谱密度针对的是功率有限的信号，例如周期信号和噪声。功率谱密度反映了单位频带内信号功率随频率的变化情况。功率谱密度与x轴包围的面积是信号的平均功率（即总能量除以信号长度$T$）。

下面列出各种量的表达式：

时域信号$\mathrm{[unit]}: x(t)=\int_{-\infty}^{+\infty}X(f)\exp(2\pi ift)\mathrm{d}f$

频域信号$\mathrm{[unit\cdot s]}: X(f)=\int_{-\infty}^{+\infty}x(t)\exp(-2\pi ift)\mathrm{d}t$

能量谱密度$\mathrm{[unit^2 \cdot s\cdot Hz^{-1}]}: |X(f)|^2$

瞬时功率$\mathrm{[unit^2 ]}: |x(t)|^2$

Parseval 定理：总能量$\mathrm{[unit^2\cdot s]} = \int_{-\infty}^{+\infty}|x(t)|^2\mathrm{d}t=\int_{-\infty}^{+\infty}|X(f)|^2\mathrm{d}f$

截取$x(t)$的一段$x_T(t)$，则 $\lim_{T \to \infty}x_T(t)=x(t)$

平均功率$\mathrm{[unit^2]}: P=\lim_{T \to \infty}\int^{+\infty}_{-\infty}|X(f)|^2\mathrm{d}f$

功率谱密度$\mathrm{[unit^2\cdot Hz^{-1}]}: S_{xx}(f)=\lim_{T \to \infty}\frac{1}{T}|X(f)|^2$

功率谱密度无法精确获得，只能采样进行谱估计。

三、功率谱密度的估计

1.周期图法

周期图法是迄今为止最经典的方法，由Arthur Schuster在1898年提出：

$$S(f)=\frac{1}{N}|\sum^{N-1}_{n=0}x_n\exp{-2\pi ift}|^2, $$

功率谱密度$$S_{xx}(f)=S(f)/\mathrm{Fs},$$ Fs为采样频率。下面的代码演示了用周期图法估算功率谱密度：

Pxx = np.abs(fft(signal))**2/N/Fs*2 #单边谱，非零频率要乘2
Pxx[0] = 0.5*Pxx[0]
freq0 = Fs*list0/N


freqs = freq0[:len(signal)//2]
Pxx = Pxx[:len(signal)//2]

plt.plot(freqs, Pxx)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power Spectral Density(unit^2/Hz)')
plt.show()

如何理解上图的纵坐标呢？

对于交流分量，信号的有效值是幅值乘以二分之根号二；对于直流分量，即频率为0，信号的有效值就是幅值。将有效值的平方乘以时长T，就约等于纵坐标的值(如果信号不加噪声，则约等于变为等于)。大家可以用上面的例子验证一下。

2.其他估计方法

周期图是信号功率谱的一个有偏估值，当信号的长度增大到无穷时，估值的方差不趋于零。Welch提出了修正的周期图功率谱密度的估计方法，它将信号分段加窗求其功率谱密度，然后做平均处理。此外，基于周期图的技术引入了一些在某些应用中不可接受的小偏差，参数方法得到了发展。参数化方法假设潜在的平稳随机过程有一个特定的结构，可以用少量的参数来描述，例如Autoregressive model (AR) 估计和Moving-average model (MA) 估计。

四、扩展

功率谱除了应用于频域，还可以用在任何的一对对偶空间上，例如宇宙学中的密度扰动功率谱是在波数空间中，温度扰动角功率谱是在球谐展开的级数l中，其中功率谱的含义，需要根据实际情况迁移领会。

这里举一个实际应用的例子。例如比较两个密度扰动场在统计上的相似程度，首先分别求两个场的功率谱和互相关功率谱，最后用“关联功率谱/根号下(功率谱1*功率谱2)”得到互相关系数。互相关系数越接近1，说明在对应波数下越相似。

Smilie Vote is loading.