2.5 `scipy.signal`：信号处理

2.5 scipy.signal：信号处理

2.5 scipy.signal: 信号处理

scipy.signal 模块是 SciPy 库中专门用于信号处理的工具箱。它提供了广泛的函数和类，用于信号的生成、滤波、分析和处理。无论你是处理音频、图像、生物信号还是其他类型的数据，scipy.signal 都能提供强大的支持。

2.5.1 信号生成

在进行信号处理之前，首先需要有信号。scipy.signal 提供了生成各种常见信号的函数。

• scipy.signal.unit_impulse(shape, idx=None, dtype=<class 'float'>): 生成单位脉冲信号（也称为狄拉克δ函数）。

  import numpy as np
  import matplotlib.pyplot as plt
  from scipy import signal
  # 生成长度为 10 的单位脉冲信号，脉冲位置在索引 5
  impulse = signal.unit_impulse(10, 5)
  plt.stem(impulse)
  plt.title("Unit Impulse Signal")
  plt.xlabel("Sample")
  plt.ylabel("Amplitude")
  plt.show()
  ​

• 其他信号生成函数: scipy.signal 还支持生成其他信号，例如正弦波、方波等。可以使用 NumPy 生成时间序列，然后结合 NumPy 的数学函数生成相应的信号。

  import numpy as np
  import matplotlib.pyplot as plt
  # 生成正弦波信号
  t = np.linspace(0, 1, 500, endpoint=False)  # 1 秒内的 500 个采样点
  f = 5  # 信号频率为 5 Hz
  sin_wave = np.sin(2 * np.pi * f * t)
  plt.plot(t, sin_wave)
  plt.title("Sine Wave Signal")
  plt.xlabel("Time (s)")
  plt.ylabel("Amplitude")
  plt.show()
  ​

2.5.2 线性时不变系统 (LTI)

scipy.signal 模块的核心概念之一是线性时不变 (LTI) 系统。 LTI 系统可以通过传递函数、状态空间方程或差分方程来描述。

• scipy.signal.TransferFunction: 用传递函数表示 LTI 系统。

  from scipy import signal
  # 定义一个传递函数：H(s) = (s+1) / (s^2 + 2s + 1)
  num = [1, 1]  # 分子系数
  den = [1, 2, 1]  # 分母系数
  system = signal.TransferFunction(num, den)
  print(system)
  ​

• scipy.signal.StateSpace: 用状态空间方程表示 LTI 系统。

  from scipy import signal
  import numpy as np
  # 定义状态空间方程
  A = [[0, 1], [-1, -2]]
  B = [[0], [1]]
  C = [[1, 0]]
  D = [[0]]
  system = signal.StateSpace(A, B, C, D)
  print(system)
  ​

• LTI 系统的分析: 可以对 LTI 系统进行各种分析，例如计算频率响应、阶跃响应等。

  import numpy as np
  import matplotlib.pyplot as plt
  from scipy import signal
  # 定义一个传递函数
  num = [1]
  den = [1, 1, 1]
  system = signal.TransferFunction(num, den)
  # 计算频率响应
  w, mag, phase = system.bode()
  plt.figure()
  plt.subplot(2, 1, 1)
  plt.semilogx(w, mag)  # 分贝刻度
  plt.title("Bode Plot - Magnitude")
  plt.ylabel("Magnitude (dB)")
  plt.subplot(2, 1, 2)
  plt.semilogx(w, phase)  # 角度刻度
  plt.title("Bode Plot - Phase")
  plt.xlabel("Frequency (rad/s)")
  plt.ylabel("Phase (degrees)")
  plt.tight_layout()
  plt.show()
  # 计算阶跃响应
  t, y = signal.step(system)
  plt.figure()
  plt.plot(t, y)
  plt.title("Step Response")
  plt.xlabel("Time (s)")
  plt.ylabel("Amplitude")
  plt.show()
  ​

2.5.3 滤波器设计

scipy.signal 提供了强大的滤波器设计工具，可以设计各种类型的滤波器，例如低通、高通、带通和带阻滤波器。

• scipy.signal.butter: 设计 Butterworth 滤波器。

  from scipy import signal
  import numpy as np
  import matplotlib.pyplot as plt
  # 设计一个 3 阶 Butterworth 低通滤波器，截止频率为 0.2 (归一化频率)
  order = 3
  cutoff_freq = 0.2
  b, a = signal.butter(order, cutoff_freq, btype='low', analog=False)
  # 频率响应
  w, h = signal.freqz(b, a, worN=8000)
  plt.figure()
  plt.plot(0.5*fs*w/np.pi, np.abs(h), 'b')
  plt.plot(cutoff_freq*fs, 0.5*np.sqrt(2), 'ro') # 截止频率点
  plt.axvline(cutoff_freq*fs, color='k') # 截止频率线
  plt.xlim(0, 0.5*fs)
  plt.title("Butterworth Lowpass Filter Frequency Response")
  plt.xlabel("Frequency (Hz)")
  plt.ylabel("Gain")
  plt.grid(True)
  plt.show()
  # 应用滤波器
  fs = 10  # 采样频率
  t = np.arange(0, 1, 1/fs)
  # 创建一个包含噪声的信号
  noise = np.random.normal(0, 0.5, len(t))
  signal_noisy = np.sin(2*np.pi*t) + noise
  # 使用滤波器
  filtered_signal = signal.lfilter(b, a, signal_noisy)
  plt.figure()
  plt.plot(t, signal_noisy, label="Noisy Signal")
  plt.plot(t, filtered_signal, label="Filtered Signal")
  plt.xlabel("Time (s)")
  plt.ylabel("Amplitude")
  plt.legend()
  plt.title("Butterworth Filter Application")
  plt.show()
  ​

• scipy.signal.firwin: 设计 FIR 滤波器。

  from scipy import signal
  import numpy as np
  import matplotlib.pyplot as plt
  # 设计一个 51 阶 FIR 低通滤波器，截止频率为 0.3 (归一化频率)
  numtaps = 51
  cutoff_freq = 0.3
  taps = signal.firwin(numtaps, cutoff_freq)
  # 频率响应
  w, h = signal.freqz(taps, worN=8000)
  plt.figure()
  plt.plot(0.5*fs*w/np.pi, np.abs(h), 'b')
  plt.axvline(cutoff_freq*fs, color='k') # 截止频率线
  plt.xlim(0, 0.5*fs)
  plt.title("FIR Lowpass Filter Frequency Response")
  plt.xlabel("Frequency (Hz)")
  plt.ylabel("Gain")
  plt.grid(True)
  plt.show()
  # 应用滤波器
  fs = 10  # 采样频率
  t = np.arange(0, 1, 1/fs)
  # 创建一个包含噪声的信号
  noise = np.random.normal(0, 0.5, len(t))
  signal_noisy = np.sin(2*np.pi*t) + noise
  # 使用滤波器
  filtered_signal = signal.lfilter(taps, 1.0, signal_noisy)
  plt.figure()
  plt.plot(t, signal_noisy, label="Noisy Signal")
  plt.plot(t, filtered_signal, label="Filtered Signal")
  plt.xlabel("Time (s)")
  plt.ylabel("Amplitude")
  plt.legend()
  plt.title("FIR Filter Application")
  plt.show()
  ​

• 滤波器类型: scipy.signal 支持各种滤波器类型，包括 Butterworth, Chebyshev (I 和 II 型), Elliptic (Cauer) 等。

2.5.4 滤波器的应用

设计好滤波器后，就可以将其应用于信号。

• scipy.signal.lfilter: 对信号应用 IIR 滤波器。

  from scipy import signal
  import numpy as np
  import matplotlib.pyplot as plt
  # 生成一个包含噪声的信号
  fs = 1000  # 采样频率
  t = np.linspace(0, 1, fs, endpoint=False)
  f = 10  # 信号频率
  noise = np.random.normal(0, 0.5, len(t))
  signal_noisy = np.sin(2 * np.pi * f * t) + noise
  # 设计一个低通 Butterworth 滤波器
  cutoff_freq = 30  # 截止频率
  order = 4  # 滤波器阶数
  b, a = signal.butter(order, cutoff_freq / (fs / 2), btype='low', analog=False)
  # 应用滤波器
  filtered_signal = signal.lfilter(b, a, signal_noisy)
  # 绘制结果
  plt.figure(figsize=(10, 6))
  plt.plot(t, signal_noisy, label='Noisy Signal')
  plt.plot(t, filtered_signal, label='Filtered Signal')
  plt.xlabel('Time (s)')
  plt.ylabel('Amplitude')
  plt.legend()
  plt.title('Filtering with lfilter')
  plt.grid(True)
  plt.show()
  ​

• scipy.signal.filtfilt: 对信号应用零相位滤波。 filtfilt 函数通过前后两次应用滤波器来消除相位失真，这在某些应用中非常重要。

  from scipy import signal
  import numpy as np
  import matplotlib.pyplot as plt
  # 生成一个包含噪声的信号
  fs = 1000  # 采样频率
  t = np.linspace(0, 1, fs, endpoint=False)
  f = 10  # 信号频率
  noise = np.random.normal(0, 0.5, len(t))
  signal_noisy = np.sin(2 * np.pi * f * t) + noise
  # 设计一个低通 Butterworth 滤波器
  cutoff_freq = 30  # 截止频率
  order = 4  # 滤波器阶数
  b, a = signal.butter(order, cutoff_freq / (fs / 2), btype='low', analog=False)
  # 应用零相位滤波
  filtered_signal = signal.filtfilt(b, a, signal_noisy)
  # 绘制结果
  plt.figure(figsize=(10, 6))
  plt.plot(t, signal_noisy, label='Noisy Signal')
  plt.plot(t, filtered_signal, label='Filtered Signal')
  plt.xlabel('Time (s)')
  plt.ylabel('Amplitude')
  plt.legend()
  plt.title('Filtering with filtfilt (Zero-Phase Filtering)')
  plt.grid(True)
  plt.show()
  ​

2.5.5 频谱分析

scipy.signal 提供了频谱分析工具，用于分析信号的频率成分。

• scipy.signal.fftconvolve: 使用 FFT 进行卷积运算。

• scipy.signal.spectrogram: 计算信号的频谱图。

  import numpy as np
  import matplotlib.pyplot as plt
  from scipy import signal
  # 生成一个信号
  fs = 1000  # 采样频率
  t = np.linspace(0, 1, fs, endpoint=False)
  f1 = 10  # 第一个频率分量
  f2 = 50  # 第二个频率分量
  signal_ = np.sin(2 * np.pi * f1 * t) + np.sin(2 * np.pi * f2 * t)
  # 计算频谱图
  f, t, Sxx = signal.spectrogram(signal_, fs)
  # 绘制频谱图
  plt.figure(figsize=(10, 6))
  plt.pcolormesh(t, f, np.log10(Sxx), shading='gouraud')
  plt.ylabel('Frequency (Hz)')
  plt.xlabel('Time (s)')
  plt.title('Spectrogram')
  plt.colorbar(label='Log Magnitude')
  plt.ylim(0, 100)  # 限制频率范围
  plt.show()
  ​

• scipy.signal.welch: 使用 Welch 方法估计功率谱密度。

  import numpy as np
  import matplotlib.pyplot as plt
  from scipy import signal
  # 生成一个信号
  fs = 1000  # 采样频率
  t = np.linspace(0, 1, fs, endpoint=False)
  f1 = 10  # 第一个频率分量
  f2 = 50  # 第二个频率分量
  noise = np.random.normal(0, 0.1, len(t))
  signal_ = np.sin(2 * np.pi * f1 * t) + np.sin(2 * np.pi * f2 * t) + noise
  # 使用 Welch 方法估计功率谱密度
  f, Pxx = signal.welch(signal_, fs, nperseg=256)
  # 绘制功率谱密度
  plt.figure(figsize=(10, 6))
  plt.semilogy(f, Pxx)
  plt.xlabel('Frequency (Hz)')
  plt.ylabel('Power Spectral Density')
  plt.title('Welch Power Spectral Density Estimate')
  plt.grid(True)
  plt.show()
  ​

2.5.6 窗函数

窗函数用于减少频谱分析中的泄漏效应。scipy.signal 提供了各种窗函数。

• scipy.signal.windows: 包含各种窗函数，例如汉宁窗、汉明窗、布莱克曼窗等。

  from scipy import signal
  import numpy as np
  import matplotlib.pyplot as plt
  # 生成一个汉宁窗
  window = signal.windows.hann(50)
  # 绘制窗函数
  plt.figure(figsize=(8, 4))
  plt.plot(window)
  plt.title("Hanning Window")
  plt.xlabel("Sample")
  plt.ylabel("Amplitude")
  plt.grid(True)
  plt.show()
  ​

2.5.7 信号卷积

• scipy.signal.convolve: 计算两个信号的线性卷积。

• scipy.signal.correlate: 计算两个信号的互相关。

  import numpy as np
  from scipy import signal
  import matplotlib.pyplot as plt
  # 定义两个信号
  signal1 = np.array([1, 2, 3, 4, 5])
  signal2 = np.array([0, 1, 0.5])
  # 计算卷积
  convolution = signal.convolve(signal1, signal2, mode='full') # full, valid, same
  print("Convolution:", convolution)
  # 计算互相关
  correlation = signal.correlate(signal1, signal2, mode='full')
  print("Correlation:", correlation)
  # 可视化
  plt.figure(figsize=(10, 5))
  plt.subplot(2, 1, 1)
  plt.plot(signal1, label='Signal 1')
  plt.plot(signal2, label='Signal 2')
  plt.title('Signals')
  plt.legend()
  plt.subplot(2, 1, 2)
  plt.plot(convolution, label='Convolution')
  plt.plot(correlation, label='Correlation')
  plt.title('Convolution and Correlation')
  plt.legend()
  plt.tight_layout()
  plt.show()
  ​

2.5.8 重采样

• scipy.signal.resample: 将信号重采样到新的采样率。

  from scipy import signal
  import numpy as np
  import matplotlib.pyplot as plt
  # 生成一个信号
  fs = 100  # 原始采样频率
  t = np.linspace(0, 1, fs, endpoint=False)
  f = 10  # 信号频率
  signal_ = np.sin(2 * np.pi * f * t)
  # 重采样到新的采样频率
  fs_new = 50  # 新的采样频率
  num = int(len(signal_) * fs_new / fs)
  signal_resampled = signal.resample(signal_, num)
  t_new = np.linspace(0, 1, len(signal_resampled), endpoint=False)
  # 绘制结果
  plt.figure(figsize=(10, 6))
  plt.plot(t, signal_, label='Original Signal')
  plt.plot(t_new, signal_resampled, label='Resampled Signal')
  plt.xlabel('Time (s)')
  plt.ylabel('Amplitude')
  plt.legend()
  plt.title('Signal Resampling')
  plt.grid(True)
  plt.show()
  ​

2.5.9 scipy.signal 模块功能总结

2.5.10 总结

scipy.signal 模块是 SciPy 库中强大的信号处理工具箱。它提供了广泛的函数和类，用于信号的生成、滤波、分析和处理。 通过学习和使用 scipy.signal，可以高效地进行各种信号处理任务。希望本文能够帮助你更好地理解和应用 scipy.signal 模块。