Python的信号处理：使用`SciPy`和`NumPy`进行信号滤波和频域分析。 - 智猿学院-前后端，数据库，人工智能，云计算等领域前沿技术讲座

Python信号处理：SciPy和NumPy的信号滤波与频域分析

大家好！今天我们来深入探讨如何使用Python中的SciPy和NumPy库进行信号处理，特别是信号滤波和频域分析。这两个库是Python科学计算的核心，为我们提供了强大的工具来处理各种信号，从音频到生物信号，再到金融数据。

1. 信号表示与NumPy

信号本质上是随时间或空间变化的物理量的表示。在数字世界中，信号被离散化和量化，最终以数字序列的形式存储。NumPy是Python中用于数值计算的基础库，它提供了高效的多维数组对象，非常适合表示和操作信号。

import numpy as np
import matplotlib.pyplot as plt

# 创建一个简单的正弦波信号
fs = 1000  # 采样频率 (Hz)
t = np.arange(0, 1, 1/fs)  # 时间向量 (0到1秒，采样频率为fs)
f = 5  # 信号频率 (Hz)
signal = np.sin(2 * np.pi * f * t)

# 绘制信号
plt.figure(figsize=(10, 4))
plt.plot(t, signal)
plt.xlabel('时间 (s)')
plt.ylabel('幅度')
plt.title('5 Hz 正弦波')
plt.grid(True)
plt.show()

在这个例子中，我们使用NumPy创建了一个时间向量t，然后使用np.sin()函数生成了一个5Hz的正弦波信号。fs代表采样频率，它决定了信号的离散化程度。采样频率必须至少是信号最高频率的两倍，以满足奈奎斯特采样定理，避免混叠现象。

2. 频域分析与快速傅里叶变换(FFT)

频域分析是将信号从时域转换到频域，从而揭示信号中包含的频率成分。快速傅里叶变换（FFT）是一种高效的算法，用于计算离散傅里叶变换（DFT）。SciPy的scipy.fft模块提供了FFT的实现。

from scipy.fft import fft, fftfreq

# 计算信号的FFT
N = len(signal)
yf = fft(signal)
xf = fftfreq(N, 1/fs)

# 绘制频谱
plt.figure(figsize=(10, 4))
plt.plot(xf, np.abs(yf)) #取绝对值，因为FFT结果是复数
plt.xlabel('频率 (Hz)')
plt.ylabel('幅度')
plt.title('信号频谱')
plt.grid(True)
plt.xlim(0, 10) #为了显示清晰，只显示0-10Hz范围内的频率
plt.show()

# 找到主频率
peak_index = np.argmax(np.abs(yf))
peak_frequency = xf[peak_index]
print(f"主频率: {peak_frequency:.2f} Hz")

这段代码首先使用fft()函数计算信号的FFT。fftfreq()函数生成与FFT输出对应的频率轴。然后，我们绘制了频谱，横轴是频率，纵轴是FFT结果的幅度。最后，我们找到幅度最大的频率分量，即主频率。注意，FFT的结果是复数，我们需要取绝对值才能得到幅度谱。

3. 信号滤波与SciPy

信号滤波是指选择性地改变信号中不同频率成分的幅度。SciPy的scipy.signal模块提供了各种滤波器设计和应用函数。

3.1 滤波器类型

常见的滤波器类型包括：

低通滤波器 (Lowpass Filter): 允许低频信号通过，衰减高频信号。
高通滤波器 (Highpass Filter): 允许高频信号通过，衰减低频信号。
带通滤波器 (Bandpass Filter): 允许特定频率范围内的信号通过，衰减其他频率的信号。
带阻滤波器 (Bandstop Filter): 衰减特定频率范围内的信号，允许其他频率的信号通过。

3.2 滤波器设计

SciPy提供了多种滤波器设计方法，例如：

巴特沃斯滤波器 (Butterworth Filter): 在通带内具有最大平坦响应。
切比雪夫滤波器 (Chebyshev Filter): 在通带或阻带内具有纹波，但具有更快的滚降。
贝塞尔滤波器 (Bessel Filter): 具有线性相位响应，适用于需要保持信号形状的应用。

from scipy.signal import butter, lfilter

def butter_lowpass(cutoff, fs, order=5):
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = butter(order, normal_cutoff, btype='low', analog=False)
    return b, a

def butter_lowpass_filter(data, cutoff, fs, order=5):
    b, a = butter_lowpass(cutoff, fs, order=order)
    y = lfilter(b, a, data)
    return y

# 创建一个包含噪声的信号
noise = 0.5 * np.random.randn(len(t))
noisy_signal = signal + noise

# 设计一个低通滤波器
cutoff_frequency = 8  # 截止频率 (Hz)
order = 5  # 滤波器阶数

# 应用低通滤波器
filtered_signal = butter_lowpass_filter(noisy_signal, cutoff_frequency, fs, order)

# 绘制结果
plt.figure(figsize=(12, 6))

plt.subplot(3, 1, 1)
plt.plot(t, noisy_signal)
plt.xlabel('时间 (s)')
plt.ylabel('幅度')
plt.title('噪声信号')
plt.grid(True)

plt.subplot(3, 1, 2)
plt.plot(t, filtered_signal)
plt.xlabel('时间 (s)')
plt.ylabel('幅度')
plt.title('滤波后的信号')
plt.grid(True)

plt.subplot(3, 1, 3)
yf_noisy = fft(noisy_signal)
xf_noisy = fftfreq(len(noisy_signal), 1/fs)
yf_filtered = fft(filtered_signal)
xf_filtered = fftfreq(len(filtered_signal), 1/fs)

plt.plot(xf_noisy, np.abs(yf_noisy), label='噪声信号')
plt.plot(xf_filtered, np.abs(yf_filtered), label='滤波后信号')
plt.xlabel('频率 (Hz)')
plt.ylabel('幅度')
plt.title('频谱比较')
plt.grid(True)
plt.xlim(0, 20)
plt.legend()

plt.tight_layout()
plt.show()

这段代码首先定义了一个butter_lowpass函数，用于设计巴特沃斯低通滤波器。 butter()函数返回滤波器的系数b和a，它们分别代表分子和分母的系数。lfilter()函数将滤波器应用于信号。我们创建了一个包含噪声的信号，然后使用低通滤波器对其进行滤波。最后，我们绘制了原始噪声信号、滤波后的信号以及它们的频谱，以便比较滤波效果。

3.3 使用 sosfilt 提高稳定性

对于高阶滤波器，直接使用 lfilter 可能会由于数值精度问题导致不稳定。sosfilt 函数使用二阶节 (Second-Order Sections) 来实现滤波，可以提高滤波器的稳定性。

from scipy.signal import butter, sosfilt, sosfreqz
import matplotlib.pyplot as plt
import numpy as np

# 生成测试信号
fs = 1000  # 采样频率
t = np.linspace(0, 1, fs, endpoint=False)  # 1秒的时间
f1, f2 = 10, 200  # 两个频率
x = 0.1 * np.sin(2*np.pi*f1*t) + np.sin(2*np.pi*f2*t) + np.random.randn(len(t)) #加入了噪声

# 设计 Butterworth 带通滤波器
lowcut = 50.0 #设置下限频率
highcut = 150.0 #设置上限频率
order = 10 # 滤波器阶数
nyq = 0.5 * fs
low = lowcut / nyq
high = highcut / nyq
sos = butter(order, [low, high], btype='band', output='sos')

# 应用滤波器
y = sosfilt(sos, x)

# 绘制结果
plt.figure(figsize=(12, 6))
plt.subplot(2, 1, 1)
plt.plot(t, x, label='原始信号')
plt.plot(t, y, label='滤波后的信号')
plt.xlabel('时间 (s)')
plt.ylabel('幅度')
plt.title('带通滤波')
plt.legend()
plt.grid(True)

# 绘制滤波器的频率响应
w, h = sosfreqz(sos, worN=8000, fs=fs)
plt.subplot(2, 1, 2)
plt.plot(w, abs(h))
plt.xlabel('频率 (Hz)')
plt.ylabel('幅度')
plt.title('滤波器频率响应')
plt.grid(True)
plt.xscale('log')
plt.tight_layout()
plt.show()

在这个例子中，我们使用 butter 函数的 output='sos' 参数来生成二阶节形式的滤波器系数。然后，我们使用 sosfilt 函数将滤波器应用于信号。 sosfreqz 用于绘制滤波器的频率响应。

4. 窗函数

在进行频域分析时，我们通常会对信号进行加窗处理，以减少频谱泄漏。频谱泄漏是指由于信号的截断，导致信号的能量扩散到其他频率上。窗函数可以平滑信号的边缘，从而减少频谱泄漏。

from scipy.signal import windows

# 创建一个汉宁窗
window = windows.hann(N)

# 将窗函数应用于信号
windowed_signal = signal * window

# 计算窗函数的FFT
yf_windowed = fft(windowed_signal)
xf_windowed = fftfreq(N, 1/fs)

# 绘制频谱
plt.figure(figsize=(10, 4))
plt.plot(xf, np.abs(yf), label='原始信号')
plt.plot(xf_windowed, np.abs(yf_windowed), label='加窗信号')
plt.xlabel('频率 (Hz)')
plt.ylabel('幅度')
plt.title('频谱比较 (加窗)')
plt.grid(True)
plt.xlim(0, 10)
plt.legend()
plt.show()

这段代码首先使用windows.hann()函数创建一个汉宁窗。然后，我们将窗函数应用于信号，并计算加窗信号的FFT。最后，我们绘制了原始信号和加窗信号的频谱，以便比较加窗的效果。常见的窗函数包括矩形窗、汉宁窗、海明窗和布莱克曼窗。不同的窗函数具有不同的特性，适用于不同的应用。

5. 信号处理流程示例：音频降噪

现在，我们来看一个更复杂的例子：使用SciPy和NumPy进行音频降噪。

import numpy as np
import scipy.io.wavfile as wavfile
from scipy.signal import butter, sosfilt

# 读取音频文件
try:
    rate, data = wavfile.read('noisy_audio.wav') #需要替换成你本地的音频文件
except FileNotFoundError:
    print("找不到音频文件。请确保'noisy_audio.wav'文件存在于当前目录下。")
    exit()

# 将音频数据转换为单声道
if len(data.shape) > 1:  # 如果是立体声
    data = np.mean(data, axis=1)

# 定义滤波器函数
def butter_bandstop(lowcut, highcut, fs, order=3):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    sos = butter(order, [low, high], btype='bandstop', output='sos')
    return sos

def bandstop_filter(data, lowcut, highcut, fs, order=3):
    sos = butter_bandstop(lowcut, highcut, fs, order=order)
    y = sosfilt(sos, data)
    return y

# 设置噪声频率范围
noise_lowcut = 500  # 噪声下限频率
noise_highcut = 1000 # 噪声上限频率

# 应用带阻滤波器
filtered_audio = bandstop_filter(data, noise_lowcut, noise_highcut, rate)

# 保存滤波后的音频
wavfile.write('filtered_audio.wav', rate, filtered_audio.astype(np.int16))  # 转换回int16类型

print("音频降噪完成！")

这个例子首先读取一个包含噪声的音频文件。然后，我们定义了一个带阻滤波器，用于衰减噪声频率范围内的信号。最后，我们将滤波器应用于音频数据，并将滤波后的音频保存到文件中。你需要将代码中的 noisy_audio.wav 替换成你本地的包含噪声的音频文件。

6. 总结与展望

今天，我们学习了如何使用SciPy和NumPy进行信号处理，包括信号的表示、频域分析和信号滤波。我们了解了FFT算法的原理，以及如何使用SciPy设计和应用各种滤波器。这些技术在音频处理、图像处理、生物信号处理等领域都有广泛的应用。掌握这些技能，能够帮助我们更好地理解和处理现实世界中的各种信号。

Python信号处理：SciPy和NumPy的信号滤波与频域分析

发表回复 取消回复

发表回复取消回复