生物医学工程 | EEG分析 - 傅里叶变换

傅立叶变换

数学定义

离散傅立叶变换定义

$F_n= \sum_{k=0}^{N-1} f_k e^{-2\pi ink/N}$

连续傅立叶变换定义

$F(v)=\int_{-\infty}^{+\infty} f(t)e^{-2\pi ivt}dt$

离散傅立叶逆变换定义

$f_k= \frac{1}{n}\sum_{n=0}^{N-1} F_n e^{2\pi ink/N}$

连续傅立叶逆变换定义

$f(t)=\int_{-\infty}^{+\infty} F(v)e^{2\pi ivt}dv$

解释和理解

傅立叶变换的目的是获得信号在频域上的分解（decomposition）。

其主要是通过构建一系列不同频率的正弦波（sine wave），并且分别与原信号的复数形式求点乘（dot product）。

信号通过傅立叶变换到频域（frequency-domain），然后通过傅立叶逆变换变换回时域（time-domain）

Q：为什么要信号的复数形式？

如果两个频率相近且相位一致的正弦波相乘，则结果的幅值最大。但如果两个波的频率一致，但相位差$\pi/2$，则直接相乘的结果为0。

我们并不希望这种情况发生，所以我们通过将信号变成复数（complex），即实数（real）部分和虚数（imaginary）部分之和。复数信号与某一个频率的正弦波点乘得到的结果不会受到相位差的影响。

通过欧拉公式（Euler’s Formula）我们可以将复数正弦波用向量（vector）表示：

$Ae^{2\pi ft}=A[\cos (2\pi ft) +i\sin(2\pi ft)]$

其中A表示幅度（即向量的长度），f表示旋转的频率

正弦波组的频率选择

根据奈奎斯特定理（Nyquist theorem）可知，我们可以测量的最大频率为(sampling frequency)/2

因此如果我们的离散傅立叶变换结果（复数）长度为N，则我们需要N个不同频率的正弦波才能完整的提取信号。

如果信号只有实数部分，则我们需要$N/2+1$个不同频率的正弦波来提取信号。

代码

fft｜快速傅立叶变换｜Fast Fourier transform

fft可以大幅度缩短傅立叶变换所需要的时间

MATLAB代码

% FFT的使用
dataFFT = fft(data); % 对data进行fft变换
% dataFFT的长度为1xN，有N/2+1个unique的频率（从0到fs/2）

% 因此可以通过下面代码绘制频域图
fVec = linspace(0,fs/2, N/2+1); % 构建频率向量
plot(fVec, abs(dataFFT(1:(N/2+1)))); % 只绘制非零的频率（positive frequency）

% 加入NFFT
dataFFT = fft(data, NFFT); % NFFT表示FFT的长度
% 如果NFFT>N，则对data末端进行zero-padding；如果NFFT<N，则将data从末端删去一定的点
fVec = linspace(0,fs/2, NFFT/2+1); % 构建频率向量
plot(fVec, abs(dataFFT(1:(NFFT/2+1)))); % 只绘制非零的频率（positive frequency）

% IFFT的使用
dataFFT = fft(data)
data = ifft(dataFFT); 

Python代码

dataFFT = np.fft.fft(data) # Perform FFT on data
fVec = np.linspace(0, fs/2, N//2+1)  # Construct frequency vector
plt.plot(fVec, np.abs(dataFFT[0:N//2+1]))  # Plot positive frequencies
plt.show()

！练习

题目：手搓一个离散傅立叶变换，信号X长度为N，采样频率（sampling frequency）为fs。

提示1：提取(N/2)+1个独立的频率，最大值为fs/2

提示2：欧拉公式：$Ae^{2\pi ft}=A[\cos (2\pi ft) +i\sin(2\pi ft)]$

可用函数：exp(t), a+bi, linspace()

N = 1000; % Length of signal
fs = 200; % Sampling frequency
time = (1:N)/fs; % Time vector
X = sin(2*pi*10*time) + 0.5*cos(2*pi*37*time); % Our signal

% Frequency information
nFrequency = N/2 + 1; % # of frequencies we can extract (including 0)
nyquist = fs/2; % Nyquist frequency
frequencies = linspace(0,nyquist,nFrequency); % Frequency vector

% Initialize Fourier coefficients
fourier = zeros(1,nFrequency);

% Fourier transform
for i=1:nFrequency
  sine_wave = exp(-1i*2*pi*frequencies(i).*time); % 构建sine族
  fourier(i) = sine_wave*X'; % 傅立叶变换为原始信号X和sine族的点乘
end

% Plot the magnitude of the Fourier coefficients
figure; plot(frequencies, abs(fourier))

import numpy as np
import matplotlib.pyplot as plt

N = 1000  # Length of signal
fs = 200  # Sampling frequency
time = np.arange(1, N+1) / fs  # Time vector
X = np.sin(2*np.pi*10*time) + 0.5*np.cos(2*np.pi*37*time)  # Our signal

# Frequency information
nFrequency = N//2 + 1  # # of frequencies we can extract (including 0)
nyquist = fs/2  # Nyquist frequency
frequencies = np.linspace(0, nyquist, nFrequency)  # Frequency vector

# Initialize Fourier coefficients
fourier = np.zeros(nFrequency, dtype=np.complex128)

# Fourier transform
for i in range(nFrequency):
    sine_wave = np.exp(-1j*2*np.pi*frequencies[i]*time)  # Construct sine wave
    fourier[i] = np.dot(sine_wave, X)  # Fourier transform as dot product of X and sine wave

# Plot the magnitude of the Fourier coefficients
plt.figure()
plt.plot(frequencies, np.abs(fourier))
plt.show()

短时傅立叶变换｜Short-time Fourier Tranform

由于傅立叶变换无法应用于不稳定的信号（unstationary signal），因此我们使用短时傅立叶变换。缩短window length可以提高空间分辨率，但是如果过小，则会导致频率分辨率太差。

Q：如果存在edge artifacts或者spectral leakage怎么办？

如果信号存在edge artifacts或者spectral leakage，则我们使用tapering对信号进行调制。常见的tapers方法有使用Hann窗、Hamming窗或者Gaussian窗。通常Hann窗可以消除edge effect。

Q：参数应该怎么设置？

time segment length：通常最好让segment能够包含2-3个cycles。segment越长 → 频率分辨率越高；segment越短 → 时间分辨率越高

overlap：通常设置为50%-90%

---

附录：

傅立叶变换理解：https://zhuanlan.zhihu.com/p/19759362

---

声明：此blog内容为上课笔记，仅为分享使用。部分图片和内容取材于课本、老师课件、网络。如果有侵权，请联系aursus.blog@gmail.com删除。