汤不热吧傅里叶变换是科学计算与信号处理领域最核心的数学工具之一。它将时域信号分解为不同频率的正弦波叠加，使我们能够在频域中分析、处理和重构信号。在Python生态中，NumPy和SciPy提供了完整而高效的傅里叶变换实现，从基本的FFT到高级的滤波器设计，覆盖了工程与科研的绝大多数需求。
本文将从傅里叶变换的基本原理出发，结合大量代码实战，系统讲解如何使用NumPy和SciPy完成频谱分析、窗口函数应用、信号滤波与滤波器设计等核心任务。所有示例均可直接运行，适合有一定Python基础的开发者快速上手。
一、傅里叶变换基础与NumPy FFT
傅里叶变换的核心思想是：任何满足一定条件的周期函数，都可以表示为不同频率、不同幅值和相位的正弦波的叠加。离散傅里叶变换（DFT）是其在计算机上的实现形式，而快速傅里叶变换（FFT）则是DFT的高效算法，将时间复杂度从O(N²)降低到O(N log N)。
NumPy的
1 | numpy.fft |
模块提供了完整的FFT实现。让我们从最基础的频谱分析开始：
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
import matplotlib.pyplot as plt
# 生成含三个频率分量的混合信号
fs = 1000 # 采样率 1000Hz
t = np.arange(0, 1, 1/fs) # 1秒时长
# 信号 = 50Hz(幅值3) + 120Hz(幅值1.5) + 300Hz(幅值0.8)
signal = 3 * np.sin(2 * np.pi * 50 * t) + \
1.5 * np.sin(2 * np.pi * 120 * t) + \
0.8 * np.sin(2 * np.pi * 300 * t)
# 执行FFT
N = len(signal)
fft_result = np.fft.fft(signal)
freq = np.fft.fftfreq(N, 1/fs)
# 取单边频谱
magnitude = np.abs(fft_result[:N//2]) / N * 2 # 乘2还原单边幅值
freq_half = freq[:N//2]
print(f"50Hz处幅值: {magnitude[np.argmin(np.abs(freq_half - 50))]:.2f}")
print(f"120Hz处幅值: {magnitude[np.argmin(np.abs(freq_half - 120))]:.2f}")
print(f"300Hz处幅值: {magnitude[np.argmin(np.abs(freq_half - 300))]:.2f}")
# 输出:
# 50Hz处幅值: 3.00
# 120Hz处幅值: 1.50
# 300Hz处幅值: 0.80
上面的代码清晰地展示了FFT的基本用法。几个关键点需要特别注意：
- 采样定理：采样率fs必须大于信号最高频率的2倍（奈奎斯特频率），否则会产生频谱混叠
- 频率分辨率：等于fs/N，即采样率除以采样点数，增加采样时长可以提高频率分辨率
- 幅值还原：FFT结果是复数，取模后需要除以N再乘2（单边谱）才能还原真实幅值
1.1 实数信号的rfft优化
对于实数信号，FFT结果具有共轭对称性——正频率和负频率的幅值相同。NumPy提供了
1 | rfft |
函数，只计算正频率部分，计算量减半：
2
3
4
5
6
7
8
9
10
fft_r = np.fft.rfft(signal)
magnitude_r = np.abs(fft_r) / N * 2
freq_r = np.fft.rfftfreq(N, 1/fs)
print(f"频率点数对比: fft={N}, rfft={len(fft_r)}")
print(f"rfft 50Hz幅值: {magnitude_r[np.argmin(np.abs(freq_r - 50))]:.2f}")
# 输出:
# 频率点数对比: fft=1000, rfft=501
# rfft 50Hz幅值: 3.00
1.2 逆变换与信号重构
使用
1 | ifft |
可以从频域数据完美还原时域信号：
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
reconstructed = np.fft.ifft(fft_result)
# 验证还原精度
max_error = np.max(np.abs(signal - reconstructed.real))
print(f"最大重构误差: {max_error:.2e}") # ~1.1e-15，浮点精度内完美还原
# 实际应用：频域滤波后重构
# 频域置零去除300Hz分量
fft_filtered = fft_result.copy()
fft_filtered[np.abs(freq - 300) < 1] = 0 # 将300Hz附近置零
fft_filtered[np.abs(freq + 300) < 1] = 0 # 负频率也要处理
filtered_signal = np.fft.ifft(fft_filtered).real
print(f"滤波后300Hz分量残余: "
f"{np.max(np.abs(filtered_signal - 3*np.sin(2*np.pi*50*t) - 1.5*np.sin(2*np.pi*120*t))):.2e}")
# 输出: ~2.2e-16，300Hz分量被成功移除
二、频谱泄漏与窗口函数
在实际应用中，我们只能截取有限长度的信号进行分析。如果截取的信号在边界处不连续（即不满足周期性），就会产生频谱泄漏——信号的能量从真实频率”泄漏”到邻近频率，导致频谱模糊。
窗口函数通过在信号两端平滑衰减到零来减少边界不连续性，从而抑制频谱泄漏。NumPy和SciPy提供了丰富的窗口函数：
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
# 生成一个非整数周期的50Hz信号（频谱泄漏场景）
fs = 256
t = np.arange(0, 1, 1/fs)
signal = np.sin(2 * np.pi * 50.5 * t) # 50.5Hz，非整数周期
# 不同窗口函数对比
windows = {
'矩形窗(无窗)': np.ones(len(t)),
'Hanning窗': np.hanning(len(t)),
'Hamming窗': np.hamming(len(t)),
'Blackman窗': np.blackman(len(t)),
'Kaiser窗(beta=8)': np.kaiser(len(t), 8)
}
print("窗口函数主瓣宽度与旁瓣衰减对比:")
print(f"{'窗口':<20} {'主瓣3dB宽度(Hz)':<20} {'第一旁瓣(dB)'}")
print("-" * 60)
for name, win in windows.items():
windowed = signal * win
fft_w = np.fft.rfft(windowed)
mag_w = 20 * np.log10(np.abs(fft_w) / np.max(np.abs(fft_w)) + 1e-10)
freq_w = np.fft.rfftfreq(len(t), 1/fs)
# 估算主瓣宽度(3dB点)
peak_idx = np.argmax(np.abs(fft_w))
peak_val = np.abs(fft_w[peak_idx])
half_power = peak_val / np.sqrt(2)
above_half = np.where(np.abs(fft_w) >= half_power)[0]
main_lobe_width = (above_half[-1] - above_half[0]) * fs / len(t)
# 估算第一旁瓣
side_lobe_region = np.abs(fft_w)[:peak_idx] if peak_idx > 5 else np.abs(fft_w)[peak_idx+5:]
if len(side_lobe_region) > 0:
first_side = np.max(side_lobe_region)
side_db = 20 * np.log10(first_side / peak_val + 1e-10)
else:
side_db = 0
print(f"{name:<20} {main_lobe_width:<20.1f} {side_db:.1f}")
窗口函数的选择需要在主瓣宽度和旁瓣衰减之间权衡：
| 窗口函数 | 主瓣宽度 | 旁瓣衰减 | 适用场景 |
|---|---|---|---|
| 矩形窗 | 最窄 | 最差(-13dB) | 频率精度优先 |
| Hanning | 中等 | 良好(-31dB) | 通用频谱分析 |
| Hamming | 中等 | 较好(-42dB) | 语音处理 |
| Blackman | 较宽 | 优秀(-58dB) | 高动态范围 |
| Kaiser | 可调 | 可调 | 灵活定制 |
2.1 短时傅里叶变换(STFT)
对于频率随时间变化的非平稳信号，全局FFT无法捕捉时变特性。短时傅里叶变换通过滑动窗口实现时频分析：
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
# 生成变频信号：前半段100Hz，后半段300Hz
fs = 2000
t = np.arange(0, 2, 1/fs)
signal = np.where(t < 1,
np.sin(2*np.pi*100*t),
np.sin(2*np.pi*300*t))
# 计算STFT
f, t_stft, Zxx = stft(signal, fs=fs, nperseg=256, noverlap=200)
print(f"STFT输出形状: {Zxx.shape}")
print(f"频率范围: {f[0]:.1f}Hz - {f[-1]:.1f}Hz")
print(f"时间范围: {t_stft[0]:.3f}s - {t_stft[-1]:.3f}s")
# 在1秒附近检查频率变化
idx_05s = np.argmin(np.abs(t_stft - 0.5))
idx_15s = np.argmin(np.abs(t_stft - 1.5))
peak_freq_early = f[np.argmax(np.abs(Zxx[:, idx_05s]))]
peak_freq_late = f[np.argmax(np.abs(Zxx[:, idx_15s]))]
print(f"0.5秒处主频: {peak_freq_early:.0f}Hz") # ~100Hz
print(f"1.5秒处主频: {peak_freq_late:.0f}Hz") # ~300Hz
三、SciPy信号滤波实战
滤波是信号处理中最常见的操作——从传感器数据中去除噪声、从音频中提取特定频段、从脑电信号中分离不同节律。SciPy的
1 | scipy.signal |
模块提供了完整的滤波工具链。
3.1 基本滤波器类型与设计
SciPy支持两大类滤波器设计方法：
- IIR（无限脉冲响应）滤波器：Butterworth、Chebyshev I/II、Elliptic等，阶数低、计算快但可能有相位失真
- FIR（有限脉冲响应）滤波器：线性相位、稳定但阶数高、计算量大
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
# 设计Butterworth低通滤波器
fs = 1000 # 采样率
cutoff = 100 # 截止频率100Hz
order = 5
b, a = butter(order, cutoff, btype='low', fs=fs)
# 查看滤波器频率响应
w, h = freqz(b, a, worN=2000, fs=fs)
response_db = 20 * np.log10(np.abs(h) + 1e-10)
cutoff_idx = np.argmin(np.abs(w - cutoff))
print(f"截止频率{cutoff}Hz处衰减: {response_db[cutoff_idx]:.1f}dB")
print(f"通带波纹(50Hz处): {response_db[np.argmin(np.abs(w-50))]:.1f}dB")
print(f"阻带衰减(200Hz处): {response_db[np.argmin(np.abs(w-200))]:.1f}dB")
# 输出:
# 截止频率100Hz处衰减: -3.0dB
# 通带波纹(50Hz处): -0.0dB
# 阻带衰减(200Hz处): -30.5dB
3.2 零相位滤波：filtfilt vs lfilter
IIR滤波器会引入相位延迟，在许多应用中这是不希望的。
1 | filtfilt |
通过前向-后向两次滤波实现零相位失真，代价是滤波器阶数等效翻倍：
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
fs = 500
t = np.arange(0, 2, 1/fs)
clean = np.sin(2*np.pi*5*t)
noise = 0.5 * np.sin(2*np.pi*50*t)
signal = clean + noise
# 设计低通滤波器（截止20Hz）
b, a = butter(4, 20, btype='low', fs=fs)
# 方法1：单向滤波（有相位延迟）
filtered_one = lfilter(b, a, signal)
# 方法2：零相位滤波（无相位延迟）
filtered_zero = filtfilt(b, a, signal)
# 对比相位偏移
mid = len(t) // 2
phase_offset_one = np.argmax(np.correlate(clean[mid-200:mid+200],
filtered_one[mid-200:mid+200], 'full')) - 400
phase_offset_zero = np.argmax(np.correlate(clean[mid-200:mid+200],
filtered_zero[mid-200:mid+200], 'full')) - 400
print(f"单向滤波相位偏移: {phase_offset_one} 采样点")
print(f"零相位滤波相位偏移: {phase_offset_zero} 采样点")
# 对比滤波效果
noise_remaining_one = np.std(filtered_one - clean)
noise_remaining_zero = np.std(filtered_zero - clean)
print(f"单向滤波残余噪声标准差: {noise_remaining_one:.4f}")
print(f"零相位滤波残余噪声标准差: {noise_remaining_zero:.4f}")
3.3 带通与带阻滤波
实际应用中，我们常常需要提取或去除特定频段。带通滤波保留指定频率范围，带阻滤波（陷波器）去除特定频率干扰：
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
fs = 360 # MIT-BIH心电数据库标准采样率
t = np.arange(0, 5, 1/fs)
# 模拟心电信号（简化）
ecu = np.sin(2*np.pi*1*t) # 基线
ecu += 0.5 * np.sin(2*np.pi*10*t) # QRS波成分
ecu += 0.3 * np.sin(2*np.pi*25*t) # 高频成分
ecu += 0.8 * np.sin(2*np.pi*60*t) # 工频干扰(60Hz)
ecu += 0.4 * np.random.randn(len(t)) # 肌电噪声
# 带通滤波：0.5-40Hz
b_bp, a_bp = butter(4, [0.5, 40], btype='band', fs=fs)
ecu_bp = filtfilt(b_bp, a_bp, ecu)
# 陷波滤波：去除60Hz工频
b_notch, a_notch = butter(4, [58, 62], btype='bandstop', fs=fs)
ecu_notch = filtfilt(b_notch, a_notch, ecu)
# 先带通再陷波
ecu_clean = filtfilt(b_notch, a_notch, ecu_bp)
print(f"原始信号SNR估算: {10*np.log10(np.var(ecu - 0.8*np.sin(2*np.pi*60*t)) / np.var(0.8*np.sin(2*np.pi*60*t))):.1f}dB")
print(f"带通滤波后60Hz幅值: {np.max(np.abs(np.fft.rfft(ecu_bp)[np.argmin(np.abs(np.fft.rfftfreq(len(t),1/fs)-60))]))/len(t)*2:.3f}")
print(f"陷波滤波后60Hz幅值: {np.max(np.abs(np.fft.rfft(ecu_notch)[np.argmin(np.abs(np.fft.rfftfreq(len(t),1/fs)-60))]))/len(t)*2:.3f}")
print(f"组合滤波后60Hz幅值: {np.max(np.abs(np.fft.rfft(ecu_clean)[np.argmin(np.abs(np.fft.rfftfreq(len(t),1/fs)-60))]))/len(t)*2:.3f}")
四、FIR滤波器设计与高级技巧
FIR滤波器具有严格的线性相位特性，在对相位敏感的应用（如音频处理、通信系统）中不可替代。SciPy提供了多种FIR设计方法：
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
fs = 1000
# 方法1：窗函数法设计低通FIR
taps_window = firwin(101, 100, fs=fs) # 101阶，截止100Hz
# 方法2：频率采样法
taps_fsamp = firwin2(101, [0, 90, 110, fs/2], [1, 1, 0, 0], fs=fs)
# 方法3：Parks-McClellan最优等波纹设计
taps_remez = remez(101, [0, 80, 120, fs/2], [1, 0], fs=fs)
# 对比过渡带宽度和阻带衰减
for name, taps in [('窗函数法', taps_window),
('频率采样法', taps_fsamp),
('等波纹法', taps_remez)]:
w, h = freqz(taps, worN=2000, fs=fs)
h_db = 20 * np.log10(np.abs(h) + 1e-10)
# 阻带衰减（150Hz以上最大值）
stopband = h_db[w > 150]
stop_atten = np.min(stopband)
# 过渡带宽度（-3dB到-40dB）
pass_edge = w[np.argmin(np.abs(h_db + 3))]
stop_edge = w[np.argmin(np.abs(h_db[:np.argmin(np.abs(h_db+40))] + 40))] if np.any(h_db < -40) else fs/2
print(f"{name}: 阻带衰减={-stop_atten:.1f}dB")
print(f"\nFIR阶数: {len(taps_window)} (线性相位，群延迟={(len(taps_window)-1)//2}采样点)")
4.1 多速率信号处理
当需要改变采样率时——如降采样减少数据量或上采样进行精细处理——SciPy提供了
1 | decimate |
、
1 | resample |
和
1 | upfirdn |
等工具：
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
# 原始信号：高采样率采集
fs_high = 10000
t = np.arange(0, 1, 1/fs_high)
signal = np.sin(2*np.pi*100*t) + 0.5*np.sin(2*np.pi*250*t)
# 降采样到1000Hz（10倍抽取，自动抗混叠滤波）
signal_dec = decimate(signal, 10, ftype='fir')
fs_dec = fs_high // 10
# 重采样到2000Hz
signal_resamp = resample(signal, int(len(signal) * 2000 / fs_high))
fs_resamp = 2000
print(f"原始: {len(signal)} 采样点, fs={fs_high}Hz")
print(f"抽取: {len(signal_dec)} 采样点, fs={fs_dec}Hz")
print(f"重采样: {len(signal_resamp)} 采样点, fs={fs_resamp}Hz")
# 验证抽取后频谱正确性
fft_orig = np.abs(np.fft.rfft(signal)) / len(signal) * 2
fft_dec = np.abs(np.fft.rfft(signal_dec)) / len(signal_dec) * 2
freq_orig = np.fft.rfftfreq(len(signal), 1/fs_high)
freq_dec = np.fft.rfftfreq(len(signal_dec), 1/fs_dec)
# 100Hz处的幅值应该保持
amp_100_orig = fft_orig[np.argmin(np.abs(freq_orig - 100))]
amp_100_dec = fft_dec[np.argmin(np.abs(freq_dec - 100))]
print(f"\n100Hz幅值 - 原始: {amp_100_orig:.3f}, 抽取后: {amp_100_dec:.3f}")
五、功率谱密度估计
功率谱密度（PSD）描述信号功率在频域的分布，是随机信号分析的核心工具。SciPy提供了周期图法和Welch方法两种主要估计方式：
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
# 生成含噪声的信号
fs = 1000
t = np.arange(0, 10, 1/fs) # 10秒长信号
np.random.seed(42)
# 有用信号 + 白噪声
signal = np.sin(2*np.pi*50*t) + 0.5*np.sin(2*np.pi*120*t)
noise = 2 * np.random.randn(len(t))
x = signal + noise
# 方法1：周期图法（方差大，不推荐直接使用）
f_pgm, Pxx_pgm = periodogram(x, fs=fs)
# 方法2：Welch方法（分段平均，方差小）
f_welch, Pxx_welch = welch(x, fs=fs, nperseg=1024, noverlap=512)
# 方法3：Welch方法使用不同窗口长度
f_welch_long, Pxx_welch_long = welch(x, fs=fs, nperseg=4096, noverlap=3072)
print("功率谱密度估计方法对比:")
print(f"{'方法':<25} {'50Hz峰值(dB)':<15} {'方差(1dB附近)'}")
print("-" * 55)
for name, f, Pxx in [('周期图法', f_pgm, Pxx_pgm),
('Welch(nperseg=1024)', f_welch, Pxx_welch),
('Welch(nperseg=4096)', f_welch_long, Pxx_welch_long)]:
# 50Hz处的PSD值
idx_50 = np.argmin(np.abs(f - 50))
peak_db = 10 * np.log10(Pxx[idx_50])
# 估算平坦区域的方差
flat_mask = (f > 200) & (f < 400)
flat_db = 10 * np.log10(Pxx[flat_mask])
variance = np.var(flat_db)
print(f"{name:<25} {peak_db:<15.1f} {variance:.2f}")
Welch方法的关键参数选择：
- nperseg：每段长度，越长频率分辨率越高但平滑效果越差
- noverlap：段间重叠，通常设为nperseg的50%-75%
- window：窗口函数，默认Hanning窗
- 工程经验：nperseg=fs*2可获得0.5Hz的频率分辨率
六、实战案例：振动信号分析与故障诊断
将前面所学技术综合应用到旋转机械振动信号分析中。这是工业界最典型的傅里叶变换应用场景之一：
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
fs = 12000 # 采样率12kHz（振动分析常用）
t = np.arange(0, 5, 1/fs)
np.random.seed(123)
# 转频25Hz（1500rpm）
rotor_freq = 25
# 正常轴承：只有转频及其谐波
normal = 2 * np.sin(2*np.pi*rotor_freq*t) # 1X
normal += 0.8 * np.sin(2*np.pi*2*rotor_freq*t) # 2X
normal += 0.3 * np.sin(2*np.pi*3*rotor_freq*t) # 3X
normal += 0.5 * np.random.randn(len(t)) # 背景噪声
# 故障轴承：出现特征频率
bpfo = 89.5 # 外圈故障频率（Ball Pass Frequency Outer）
bpfi = 130.5 # 内圈故障频率
fault = normal.copy()
fault += 1.5 * np.sin(2*np.pi*bpfo*t) # 外圈故障特征
fault += 0.8 * np.sin(2*np.pi*bpfi*t) # 内圈故障特征
# 步骤1：带通滤波去除低频漂移和高频噪声
b, a = butter(5, [10, 1000], btype='band', fs=fs)
normal_filt = filtfilt(b, a, normal)
fault_filt = filtfilt(b, a, fault)
# 步骤2：Welch功率谱估计
f_n, Pxx_n = welch(normal_filt, fs=fs, nperseg=4096)
f_f, Pxx_f = welch(fault_filt, fs=fs, nperseg=4096)
# 步骤3：特征提取
print("===== 振动分析诊断报告 =====")
print(f"转频: {rotor_freq}Hz")
print(f"外圈故障频率(BPFO): {bpfo}Hz")
print(f"内圈故障频率(BPFI): {bpfi}Hz")
print()
# 检测特征频率处的幅值
for name, f_arr, Pxx_arr in [('正常', f_n, Pxx_n), ('故障', f_f, Pxx_f)]:
amp_1x = 10*np.log10(Pxx_arr[np.argmin(np.abs(f_arr - rotor_freq))])
amp_bpfo = 10*np.log10(Pxx_arr[np.argmin(np.abs(f_arr - bpfo))])
amp_bpfi = 10*np.log10(Pxx_arr[np.argmin(np.abs(f_arr - bpfi))])
print(f"{name}轴承:")
print(f" 1X({rotor_freq}Hz): {amp_1x:.1f}dB")
print(f" BPFO({bpfo}Hz): {amp_bpfo:.1f}dB")
print(f" BPFI({bpfi}Hz): {amp_bpfi:.1f}dB")
# 步骤4：包络分析（用于冲击特征检测）
from scipy.signal import hilbert
analytic = hilbert(fault_filt)
envelope = np.abs(analytic)
# 包络谱
f_env, Pxx_env = welch(envelope, fs=fs, nperseg=4096)
print(f"\n包络谱BPFO处峰值: {10*np.log10(Pxx_env[np.argmin(np.abs(f_env - bpfo))]):.1f}dB")
print("(包络谱能有效凸显周期性冲击特征)")
七、性能优化与大规模数据处理
处理大规模信号数据时，FFT的计算效率至关重要。以下是几个实用的性能优化技巧：
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
# 1. 利用FFT对2的幂次长度的优化
sizes = [99991, 100000, 131072] # 质数、普通数、2的幂
for n in sizes:
data = np.random.randn(n)
start = time.perf_counter()
for _ in range(100):
np.fft.fft(data)
elapsed = time.perf_counter() - start
print(f"N={n:<8} 100次FFT耗时: {elapsed:.3f}s")
# 2. 使用rfft处理实数信号
data = np.random.randn(1000000)
start = time.perf_counter()
_ = np.fft.fft(data)
t_fft = time.perf_counter() - start
start = time.perf_counter()
_ = np.fft.rfft(data)
t_rfft = time.perf_counter() - start
print(f"\n100万点: fft={t_fft:.4f}s, rfft={t_rfft:.4f}s, 加速比={t_fft/t_rfft:.2f}")
# 3. 分块处理超长信号
chunk_size = 2**20 # 1M点一块
total_samples = 10 * chunk_size # 10M点
data = np.random.randn(total_samples)
start = time.perf_counter()
results = []
for i in range(0, total_samples, chunk_size):
chunk = data[i:i+chunk_size]
results.append(np.fft.rfft(chunk))
spectrum = np.mean(np.abs(results), axis=0)
elapsed = time.perf_counter() - start
print(f"\n分块处理10M点: {elapsed:.3f}s")
# 4. 预分配与就地操作
output = np.empty(500001, dtype=complex)
start = time.perf_counter()
np.fft.rfft(data[:1000000], out=output)
print(f"就地FFT: {time.perf_counter()-start:.4f}s")
性能优化的核心原则总结：
- 优先使用
1rfft
处理实数信号，计算量减半
- 数据长度尽量凑成2的幂次，可提升2-5倍速度
- 超长信号采用分块+平均策略，兼顾速度和内存
- 使用
1scipy.fft
替代
1numpy.fft可自动利用多核并行
- 实时场景中考虑
1pyFFTW
库，它使用FFTW算法并支持wisdom预规划
八、总结与最佳实践
本文系统介绍了Python科学计算中傅里叶变换与信号处理的完整工作流。以下是核心要点回顾：
- 频谱分析：使用
1np.fft.rfft
处理实数信号，注意幅值还原（除以N乘2）和频率分辨率（fs/N）
- 频谱泄漏：非整数周期截断导致泄漏，选择合适窗口函数（Hanning通用、Blackman高动态范围、Kaiser灵活可调）
- 时频分析：STFT通过
1scipy.signal.stft
实现，适用于非平稳信号
- 滤波设计：IIR用
1butter
等函数设计，FIR用
1firwin/
1remez设计，零相位滤波优先选
1filtfilt - 功率谱估计：Welch方法是工业标准，合理选择nperseg平衡分辨率与方差
- 性能优化：2的幂次长度、rfft、分块处理、scipy.fft并行
在实际工程中，建议遵循以下工作流：先对原始信号做目视检查和统计描述 → 应用合适的窗口函数 → FFT频谱分析定位特征频率 → 根据频域信息设计滤波器 → 滤波后再次频谱验证 → 如有需要进入时频分析或包络分析。这种从全局到局部、从时域到频域的系统分析方法，能够有效解决绝大多数信号处理问题。