傅里叶变换和功率谱密度

\(\)

一 、离散傅里叶变换

假设有信号\(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中,其中功率谱的含义,需要根据实际情况迁移领会。

 

Smilie Vote is loading.
✿ 阅读数:842  分类:文章

发表回复

您的电子邮箱地址不会被公开。 必填项已用*标注

Captcha Code