Scipy 是基于 Numpy 的科学计算库，用于数学、科学、工程学等领域

包含的模块有最优化、线性代数、积分、插值、特殊函数、快速傅里叶变换、信号处理和图像处理、常微分方程求解和其他科学与工程中常用的计算

信号处理

1
2
import numpy as np
import matplotlib.pyplot as plt

正弦波表达：
$$
x(t) = A_1\sin (2\pi f_1t+\phi_1)+A_2\sin (2\pi f_2t+\phi_2)
$$
添加高斯白噪声：noise_signal = signal + noise, 其中noise$\sim N(0,\sigma^2) $

快速傅里叶变换

1
from scipy import fft

傅里叶变换将时域信号变换到频域，分析不同频率分量

计算机处理的是离散采样信号，DFT 是傅里叶变换的离散形式

FFT是DFT的高效算法，将计算复杂度从$O(N^2)$降至$O(N\log N)$，极大提升计算速度

输出结果为复数数组，幅度：频率成分的强度(np.abs(x)) 相位：频率成分的偏移(np.angle(X))

numpy也有fft模块，两者的接口相同，但是只能用CPU，不支持并行或GPU，尽量用scipy库内的

信号参数：采样率$f_s$ 信号频率$f$ 采样时间$T$ 采样点数$N=f_s\cdot T$

频率分辨率：$\Delta f= f_s/N$，当信号频率不是$\Delta f$的整数倍时，能量会“泄露”到相邻频段，详见频谱泄露

奈奎斯特采样频率：采样频率至少是信号最高频率的两倍，才能无失真地重建原始信号

对于实数输入信号，频域结果具有共轭对称性：

  • 后一半频率是前一半的镜像(忽略 Nyquist 频率点)
  • 实际分析时通常取前半部分(单侧频谱)

可以只计算正频率部分，输出长度减半：

1
2
3
4
5
6
7
# 直接保留正部分
fft.rfft(sig)    
fft.rfftfreq(n,d) 
# 普通写法
fft.fft(sig)[:N//2]
fft.fftfreq(n,d)[:N//2]
# 通过`fftshift`可以换成[-fs/2...fs/2]的对称频率谱

$n$为信号长度，$d$为采样间隔，返回对应傅里叶变换的频率轴

默认的FFT函数只做了求和，所以计算出来的幅值会比实际大N倍

如果只取正频率部分，需要乘上2/N做归一化，归一化是为了保证数值物理意义正确

1
magnitude = 2/N * np.abs(signal_fft)

一维傅里叶变换：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
f1, f2, fs, T= 50, 200, 500, 1.0
N = int(fs*T) # 采样点数
t = np.linspace(0, T, N, endpoint=False)
signal1 = 0.7 * np.sin(2 * np.pi * f1 * t) + 0.5 * np.sin(2 * np.pi * f2 * t)
# 计算傅里叶变换
signal_fft = fft.rfft(signal1)
# 获取频率轴
freqs = fft.rfftfreq(N, 1/fs)
# 计算归一化幅度谱
magnitude = 2/N * np.abs(signal_fft)
reconstructed_signal = fft.ifft(signal_fft)
plt.subplot(2, 1, 1)
plt.plot(t[:100], signal1[:100], 'r')
plt.title("Time Domain Signal")
plt.xlabel("Time [s]")
plt.xlabel("Time [s]")
plt.ylabel("Amplitude")

plt.subplot(2, 1, 2)
plt.plot(freqs, magnitude, 'g')
plt.title("Magnitude Spectrum")
plt.xlabel("Frequency [Hz]")
plt.ylabel("Amplitude")
plt.tight_layout()
plt.show()

二维傅里叶：

在图像处理分析中一般不选择归一化，因为主要目的是「看结构而不是看数值」

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
# 创建一个二维信号(例如图像)
image = np.random.random((256, 256))  # 256x256的随机图像
# 对二维信号进行傅里叶变换
fft_image = fft.fft2(image)  # 不用rfft2，无法同时双轴处理
# 移动频率成分，使得零频率部分位于中心
fft_image_shifted = fft.fftshift(fft_image)
# 计算幅度谱(Magnitude Spectrum)
magnitude_spectrum = np.abs(fft_image_shifted) # 未归一化
# 计算幅度谱的对数，便于显示
log_magnitude_spectrum = np.log(magnitude_spectrum + 1)
# 可视化频谱
plt.figure(figsize=(6, 6))
plt.imshow(log_magnitude_spectrum, cmap='gray')
plt.title("Magnitude Spectrum (Log Scale)")
plt.colorbar()
plt.show()

归一化场景分析：

场景 是否需要归一化 原因 / 目标
频谱图展示(图像/信号) 不需要 只看频率分布，不关心绝对幅度
未归一化+log 更直观
信号幅值恢复 需要 FFT默认放大N倍，必须除以N才能得到真实幅度
功率/能量分析 需要 保证时域能量与频域能量一致
滤波器设计/频域运算 需要 滤波器幅值必须落在理论范围内，否则滤波结果能量失真
逆变换恢复图像/信号 看情况 如果正向和逆向保持一致就不必管；如果只做正向归一化，需要注意逆变换对应的系数

滤波器设计

1
from scipy import signal

以巴特沃斯滤波器的函数举例，各个滤波器之间的参数设置基本一致

1
butter(N, Wn, btype='low', analog=False, output='ba', fs=None)
  • N：滤波器阶数，阶数越高，滤波越陡峭
  • Wn：截止频率，是一个长度为2的序列，单个或区间上下限
  • btype：滤波器类型，包括低通滤波器(low)、高通滤波器(high)、带通滤波器(band)、带阻滤波器(bandstop)，默认为low
  • analog：默认返回一个数字滤波器，如果为True则返回一个模拟滤波器
  • output：输出，默认'ba'，可选zpk和sos
  • fs：系统的采样频率，如果不提供会报错

也可以写成：

1
butter(N, Wn/fs, btype='low', analog=False, output='ba')

IIR滤波器：

巴特沃斯滤波器：常见的滤波器，它的特点是保持平滑的频率响应，没有波动，适用于那些要求频率响应平稳的应用

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
fs = 500  # 采样频率 500Hz
lowcut = 100.0  # 截止频率 100Hz

# 使用巴特沃斯滤波器设计一个4阶低通滤波器
b, a = signal.butter(4, lowcut, btype='low', fs=fs)
# 打印滤波器系数
print(f"滤波器系数 b: {b}") # 滤波器的分子系数
print(f"滤波器系数 a: {a}") # 滤波器的分母系数
# 计算频率响应
w, h = signal.freqz(b, a, worN=2000, fs=fs)
# 绘制滤波器频率响应
plt.plot(w, abs(h), 'b')
plt.title('巴特沃斯低通滤波器频率响应')
plt.xlabel('频率 [Hz]')
plt.ylabel('幅度')
plt.grid()
plt.show()

切比雪夫滤波器：能够提供更陡峭的频率响应，但它会引入波动(波纹)，可以用于需要锐利频率选择的场景

1
2
3
4
5
6
7
8
9
10
11
12
13
rp = 1  # 通带最大波纹(dB)
lowcut = 100.0  # 截止频率 100Hz
# 使用切比雪夫I型滤波器设计一个4阶低通滤波器
b, a = signal.cheby1(4, rp, lowcut, fs=fs)
# 计算频率响应
w, h = signal.freqz(b, a, worN=2000, fs=fs)
# 绘制滤波器频率响应
plt.plot(w, abs(h), 'r')
plt.title('切比雪夫I型低通滤波器频率响应')
plt.xlabel('频率 [Hz]')
plt.ylabel('幅度')
plt.grid()
plt.show()

二阶陷波和谐振滤波器：

1
2
3
4
5
6
7
8
9
10
11
fs, f0, Q = 500, 100, 30 # f0为陷波/谐振频率
b_notch, a_notch = signal.iirnotch(f0, Q, fs=fs)   # 直接传f0,需要品质因子
b_peak, a_peak = signal.iirpeak(f0, Q, fs=fs)
# 计算频率响应
w1, h1 = signal.freqz(b_notch, a_notch, worN=2000, fs=fs)
w2, h2 = signal.freqz(b_peak, a_peak, worN=2000, fs=fs)
# 绘制滤波器频率响应
plt.subplot(2, 1, 1)
plt.plot(w1, abs(h1), 'r')
plt.subplot(2, 1, 2)
plt.plot(w2, abs(h2), 'g')

FIR 滤波器：

可以通过选择不同的窗函数来实现低通、高通等滤波效果

1
2
3
4
5
6
7
8
9
10
11
12
# 使用窗函数法设计一个低通FIR滤波器
numtaps = 101  # 滤波器的长度(系数个数)
b = signal.firwin(numtaps, lowcut, pass_zero=True, fs=fs)
# 计算频率响应
w, h = signal.freqz(b, worN=2000, fs=fs)
# 绘制滤波器频率响应
plt.plot(w, abs(h), 'g')
plt.title('FIR低通滤波器频率响应')
plt.xlabel('频率 [Hz]')
plt.ylabel('幅度')
plt.grid()
plt.show()

滤波操作

常见的两个函数：

  • scipy.signal.lfilter：使用递归公式对输入信号进行滤波

  • scipy.signal.filtfilt：对信号进行零相位滤波，即滤波操作前后都进行一次

1
2
filtered = signal.lfilter(b, a, noise_signal)
filtered = signal.filtfilt(b, a, noise_signal)
特性 lfilter filtfilt
滤波方向 单向滤波(从前到后) 双向滤波(从前到后 + 从后到前)
相位延迟 引入相位延迟(有延迟) 无相位延迟(零相位)
输出信号 会有信号的末尾和开始部分的失真 保持原始信号的时域特性，不会有失真
计算性能 计算速度较快，适用于实时信号处理 需要更多计算，适用于离线处理，较慢
使用场景 实时滤波，延迟容忍的应用 高精度滤波，要求保持信号时域特性 
1
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
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
"""
生成信号(1 kHz + 2 Khz + 50 Hz + 噪声)，进行FFT分析，
应用带阻(50 Hz)+ 巴特沃斯低通(1.5kHz)滤波，
绘制并保存结果
"""
import numpy as np
from pathlib import Path
import matplotlib.pyplot as plt
import scipy.signal as signal
import scipy.fft as fft

rng = np.random.default_rng(42)
OUTPUT_DIR = "res"

def TimeAndFFT_plot(t_seg, sig, fs, title):
    N = len(sig)
    seg_len = len(t_seg)
    freq = fft.rfftfreq(N, 1/fs)
    magnitude = 2/N*np.abs(fft.rfft(sig))
    fig, (ax1,ax2) = plt.subplots(2,1)
    ax1.plot(t_seg, sig[:seg_len], 'r--')
    ax1.set_title(title+" Time-domain")
    ax1.set_xlabel("Time (s)")
    ax1.set_ylabel("Amplitude")
    ax1.grid()

    ax2.plot(freq, magnitude, 'g')
    ax2.set_title(title+" FFT")
    ax2.set_xlabel("Frequency (Hz)")
    ax2.set_ylabel("Magnitude")
    ax2.set_xlim(0,2200)
    ax2.grid()
    plt.tight_layout()
    savepath = Path(OUTPUT_DIR)/f"{title}.png"
    plt.savefig(savepath, dpi=300)
    plt.close()

if __name__ == '__main__':
    # ==== 参数设置 ====
    Path(OUTPUT_DIR).mkdir(parents=True, exist_ok=True)
    fs = 10000
    duration = 1
    N = int(fs*duration)
    t = np.arange(N)/fs
    # ==== 信号组成 ====
    f = 1000
    f_inter = 50
    amp1 = 0.7
    amp2 = 0.6
    amp_inter = 0.5
    rand_std = 0.5

    # ==== 信号产生 ====
    signal_clean = amp1*np.sin(2*np.pi*f*t) + amp2*np.sin(2*np.pi*2*f*t)+amp_inter*np.sin(2*np.pi*f_inter*t)
    noise = 0.2*rng.normal(0, rand_std, N)
    signal_noise = signal_clean + noise

    # 用于时域可视化的短时段（前0.01秒）
    seg_len = int(0.01*fs)
    t_seg = t[:seg_len]

    # ==== 绘制时域图和频谱图 ====
    TimeAndFFT_plot(t_seg, signal_clean, fs, "original")
    TimeAndFFT_plot(t_seg, signal_noise, fs, "noisy")

    # ==== 滤波器设计 ====
    lowcut = 1500
    b_lp, a_lp = signal.butter(4, lowcut, btype='low', fs=fs) # 低通滤波器
    f0, Q = 50.0, 30
    w0 = f0/(fs/2)
    b_notch, a_notch = signal.iirnotch(w0,Q)

    # ==== 滤波操作 ====
    sig_lp = signal.filtfilt(b_lp, a_lp, signal_noise)
    sig_notch = signal.filtfilt(b_notch, a_notch, signal_noise)
    sig_lp_notch = signal.filtfilt(b_notch, a_notch, sig_lp)

    # ==== 绘制时域图和频谱图 ====
    TimeAndFFT_plot(t_seg, sig_lp, fs, "lowpass")
    TimeAndFFT_plot(t_seg, sig_notch, fs, "notch")
    TimeAndFFT_plot(t_seg, sig_lp_notch, fs, "lowpass+notch")

卷积

数学公式：
$$
y[n]=(x*h)[n]=\sum_{k=-\infty}^\infty x[k]\cdot h[n-k]
$$
信号$x[n]$和滤波器$h[n]$的“滑动加权叠加”

移动平均滤波器(最简单的滤波器，等价于低通)：
$$
y[n]=\frac{1}{M}\sum_{k=0}^{M-1}x[n-k]
$$
本质上就是卷积核 $h[n] = [1/M, 1/M, …, 1/M]$

虽然numpy库内也有卷积函数，但是大部分情况下还是选用scipy库的卷积函数，对比如下：

特性 np.convolve scipy.signal.convolve
支持维度 仅限一维 (1D) 支持任意维度(1D/2D/3D)
计算速度 小型一维数据较快 大型或多维数据优化更好
功能扩展 基础功能 支持复数、FFT加速、边界处理等高级选项 
1
signal.convolve(in1, in2, mode='same', method='auto')

mode(输出尺寸)：

  • 'full'：完整线性卷积[默认]
  • 'same'：输出与in1同形，且在full结果中居中对齐
  • 'valid'：输出仅包含那些不依赖于零填充的元素(要求in1或in2至少一个在每个维度上必须至少与另一个一样大)

method(计算方法)：

  • 'direct'：卷积直接由和来确定，即卷积的定义(适合 kernel 很短或只需少量输出点时)
  • 'fft'：通过 FFT 计算(内部调用fftconvolve)
  • 'auto'：根据对哪种方法更快进行的估计自动选择直接方法或傅里叶方法[默认]
1
2
h = np.ones(M)/M # M越大平滑效果越好
smoothed = signal.convolve(signal1, h, 'same')

频谱泄露

频谱泄漏是指在进行傅里叶变换时，由于信号截断或周期化造成的频谱畸变现象，是离散傅里叶变换DFT中一个非常常见的现象

简单来说就是原本应该集中在一个频率点上的能量，由于“泄漏”到了其他频率点上，导致频谱变得模糊不清

  • 理想频谱：单一频率的冲击函数
  • 泄露后的频谱：主瓣拓宽 + 旁瓣拖尾

为什么会发生频谱泄漏？

DFT假设信号是周期性的，且采样窗口是信号的精确周期

如果信号的频率刚好是采样频率分辨率的整数倍，那么周期延拓时首尾能够平滑拼接，频谱会集中在对应的频率点上

如果不是整数倍，则首尾无法拼接，周期性假设失效，等效于信号被突然截断，相当于乘了一个矩形窗，而在频域上，乘法相当于卷积

矩形窗的频谱是一个sinc函数，主瓣窄但旁瓣高，在频域里就会产生旁瓣扩展，能量扩散到多个频率点上，这就是频谱泄漏

本质：有限截断信号等效于矩形窗乘积 → 频域产生卷积 → 能量扩散

举个例子：

$f_s=1000Hz,N=1000\rightarrow \Delta f=f_s/N=1Hz$

如果信号频率为50Hz，对齐FFT频点，那么频谱主要集中在50Hz处；

如果信号频率为50.5Hz，不在整数频点，那么频谱主瓣在50.5Hz处，但能量会扩散到周围很多频率点，看起来拖尾

解决方法：

  • 增加采样点数N：提高频率分辨率$\Delta f= f_s/N$，让频率更接近FFT栅格
  • 加窗函数：用汉宁窗、汉明窗、布莱克曼窗等替代矩形窗，减少旁瓣泄漏，但会增宽主瓣
  • 零填充：在进行FFT变换之前，在信号末尾添加若干个零，可以提高频谱显示的平滑度(可以提高频谱分辨率，但并不能消除频谱泄漏)

实际应用要在“泄漏小”和“分辨率高”之间折中

当采样率不变的情况下，提高频率分辨率的方法就是增加采样时间

窗函数

窗函数通过在信号两端平滑过渡到零，减少边界不连续性，从而：

  1. 抑制频谱泄露(降低旁瓣)
  2. 提高频率分辨率(窄化主瓣)

常用评价指标：

  • 主瓣宽度：影响能否分辨紧邻的频率
  • 旁瓣电平：决定泄漏到远处的能量大小(越低越好)
窗函数 主瓣宽度 旁瓣峰值(dB) 常用场景
矩形窗 最窄 -13 瞬时分析
Hann窗(汉宁窗) 中等 -31 通用场景
Hamming窗(海明窗) 稍窄 -41 音频处理
Blackman窗(布莱克曼窗) 宽 -58 高精度测量 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
import numpy as np
from scipy.signal import windows

N = 1024  # 窗长度
# 生成窗函数
hanning = windows.hann(N)        # 汉宁窗
hamming = windows.hamming(N)     # 海明窗
blackman = windows.blackman(N)   # 布莱克曼窗
rect = np.ones(N)                # 矩形窗

# 可视化窗函数
plt.figure()
plt.plot(hanning, label='Hann')
plt.plot(hamming, label='Hamming')
plt.plot(blackman, label='Blackman')
plt.plot(rect, '--', label='Rectangular')
plt.legend(); plt.title('常用窗函数'); plt.grid(True)

案例

一个50Hz正弦波，采样率1000，采样215点(周期不再整数倍，出现泄露)

对比不同窗函数的频谱，并采用零填充使得频谱平滑

1
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
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
from scipy.signal import windows
import numpy as np
import matplotlib.pyplot as plt
# 设置中文字体支持
plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei', 'KaiTi', 'SimSun']
plt.rcParams['axes.unicode_minus'] = False
palette = plt.get_cmap('Set1')

# ==== 参数设置 ====
fs = 1000   # 采样率 Hz
f = 50      # 信号频率 Hz
N = 215     # 采样点数(故意不匹配) 分辨率为1/T=4.65Hz
fft_size = 1024 # 零填充，插值让曲线更平滑，建议 4~8 倍且为2的幂次
t = np.arange(N) / fs

# ==== 生成信号 ====
signal = np.sin(2 * np.pi * f * t)

# ==== 窗函数字典 ====
windows_dict = {
    "矩形窗": np.ones(N),
    "汉宁窗": windows.hann(N),
    "海明窗": windows.hamming(N),
    "布莱克曼窗": windows.blackman(N)
}

# ==== FFT 计算函数 ====
def compute_fft(sig, win, fs):
    sig_win = sig * win
    # 会自动在信号后面补零到fft_size点，再做 FFT，实现零填充
    # 使得输出曲线更平滑，但分辨率不变
    fft_vals = np.fft.rfft(sig_win, fft_size)
    freqs = np.fft.rfftfreq(fft_size, 1/fs)
    # 修正幅度：除以窗口均值，保证幅度不会因加窗衰减
    magnitude = np.abs(fft_vals) / (N * np.mean(win)) * 2
    return freqs, magnitude, sig_win

# ==== 存储结果 ====
results = {}
for name, win in windows_dict.items():
    freqs, mag, sig_win = compute_fft(signal, win, fs)
    results[name] = {"freqs": freqs, "mag": mag, "sig": sig_win}

# ==== 绘制时域信号 ====
plt.figure()
plt.title("不同窗函数的时域信号")
for name, res in results.items():
    plt.plot(t[N//2-50:N//2+50], res["sig"][N//2-50:N//2+50], label=name)
plt.xlabel("时间 (s)")
plt.ylabel("幅度")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

# ==== 分图绘制频谱对比(对数 dB) ====
plt.figure()
# enumerate(..., 1) 给循环编号
for i, (name, res) in enumerate(results.items(), 1):
    # 先归一化，让主瓣的最大值为1，再转到dB(20*log10(...))
    mag_db = 20 * np.log10(res["mag"] / np.max(res["mag"]) + 1e-12)
    plt.subplot(2, 2, i)
    plt.plot(res["freqs"], mag_db, color=palette(i))
    plt.title(f"{name} 频谱 (dB)")
    plt.xlim(0, 200)
    plt.ylim(-120, 5)
    plt.xlabel("频率 (Hz)")
    plt.ylabel("相对幅度 (dB)")
    plt.grid(True)

plt.tight_layout()
plt.show()

# ==== 绘制频域对比(对数 dB)====
plt.figure()
plt.title("不同窗函数的频谱泄漏对比")
for name, res in results.items():
    mag_db = 20 * np.log10(res["mag"] / np.max(res["mag"]) + 1e-12)
    plt.plot(res["freqs"], mag_db, label=name)
plt.axvline(x=f, color='r', linestyle='--', alpha=0.5, label=f"理论频率 {f}Hz")
plt.xlim(20, 80)
plt.ylim(-80, 5)
plt.xlabel("频率 (Hz)")
plt.ylabel("相对幅度 (dB)")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

音频处理

1
from scipy.io import wavfile

.wav 文件是 PCM 编码的无压缩音频格式

主要信息：

  • 采样率(Hz)：每秒采样次数(CD 音质是 44100 Hz)
  • 位深：每个采样点的精度(16 位、32 位)
  • 声道数：单声道 / 双声道

读入音频：

1
2
3
fs, data = wavfile.read("audio.wav")
if data.ndim > 1:
    data = data[:,0]  # 单声道

写出音频：

写 16-bit PCM(常见、兼容)，要先归一化并转换

1
2
3
4
5
6
def int16_data(data):
    data_normalized = data/ np.max(np.abs(data))
    data_int16 = (data_normalized * 32767).astype(np.int16)
    return data_int16
pcm16 = int16_data(data)
wavfile.write('sine_int16.wav', fs, pcm16)