信号处理与控制系统完全指南：从傅里叶变换到 PID 控制

• 1. 信号与系统基础
  • 1.1 信号的分类
  • 1.2 基本信号
  • 1.3 线性时不变(LTI)系统
  • 1.4 卷积
• 2. 傅里叶级数 (Fourier Series)
  • 2.1 周期信号的傅里叶级数
  • 2.2 方波(Square Wave)示例
  • 2.3 帕塞瓦尔定理
• 3. 傅里叶变换 (Fourier Transform)
  • 3.1 连续傅里叶变换 (CTFT)
  • 3.2 重要变换对
  • 3.3 傅里叶变换的性质
  • 3.4 用 Python 计算傅里叶变换
• 4. 离散傅里叶变换 (DFT) & FFT
  • 4.1 离散时间傅里叶变换 (DTFT)
  • 4.2 DFT 的定义
  • 4.3 FFT 算法
  • 4.4 频谱分析示例
  • 4.5 采样定理 (Nyquist-Shannon)
• 5. 拉普拉斯变换 (Laplace Transform)
  • 5.1 定义
  • 5.2 重要变换对
  • 5.3 传递函数 (Transfer Function)
  • 5.4 二阶系统标准型
• 6. Z 变换
  • 6.1 定义
  • 6.2 Z 变换与拉普拉斯变换的关系
  • 6.3 稳定性分析
  • 6.4 重要 Z 变换对
• 7. 数字滤波器设计
  • 7.1 FIR vs IIR 滤波器
  • 7.2 滤波器类型
  • 7.3 滤波器设计函数
  • 7.4 用 scipy.signal 设计滤波器
  • 7.5 FIR 滤波器设计 (窗函数法)
• 8. 控制系统基础
  • 8.1 开环 vs 闭环控制
  • 8.2 稳定性: 劳斯-赫尔维茨判据
  • 8.3 波特图 (Bode Plot)
• 9. PID 控制器
  • 9.1 PID 控制动作
  • 9.2 齐格勒-尼科尔斯调参法
  • 9.3 用 Python 模拟 PID
• 10. 状态空间表示 (State Space)
  • 10.1 状态方程
  • 10.2 可观测性与可控性
  • 10.3 线性二次调节器 (LQR)
• 11. AI/ML 与信号处理
  • 11.1 用于一维信号的 CNN
  • 11.2 STFT 与 Mel 频谱图
  • 11.3 信号处理 + 深度学习的应用案例
• 12. 小测验
• 参考文献

1. 信号与系统基础

1.1 信号的分类

信号(Signal)是传递信息的物理量的变化。基本分类如下。

连续时间信号 vs 离散时间信号

连续时间信号(CT Signal)在所有时间 $t \in \mathbb{R}$ 上都有定义。例：模拟语音信号 $x(t) = A\cos(\omega_0 t + \phi)$

离散时间信号(DT Signal)仅在整数索引 $n \in \mathbb{Z}$ 上有定义。例：数字音频采样 $x[n]$

能量信号与功率信号

信号的能量 $E$ 与平均功率 $P$:

$E = \int_{-\infty}^{\infty} |x(t)|^2 \, dt$

$P = \lim_{T \to \infty} \frac{1}{2T} \int_{-T}^{T} |x(t)|^2 \, dt$

若 $E < \infty$ 则为能量信号，若 $0 < P < \infty$ 则为功率信号。

1.2 基本信号

单位冲激(Unit Impulse)

$\delta(t) = \begin{cases} \infty & t = 0 \\ 0 & t \neq 0 \end{cases}, \quad \int_{-\infty}^{\infty} \delta(t) \, dt = 1$

筛选(Sifting)性质：$\int_{-\infty}^{\infty} x(t)\delta(t - t_0) \, dt = x(t_0)$

单位阶跃(Unit Step)

$u(t) = \begin{cases} 1 & t \geq 0 \\ 0 & t < 0 \end{cases}$

成立关系式 $u(t) = \int_{-\infty}^{t} \delta(\tau) \, d\tau$。

1.3 线性时不变(LTI)系统

LTI 系统满足以下两条核心性质。

• 线性(Linearity): $\mathcal{T}\{ax_1(t) + bx_2(t)\} = a\mathcal{T}\{x_1(t)\} + b\mathcal{T}\{x_2(t)\}$
• 时不变性(Time Invariance): $y(t) = \mathcal{T}\{x(t)\} \Rightarrow \mathcal{T}\{x(t-t_0)\} = y(t-t_0)$

1.4 卷积

LTI 系统的输出可以表示为输入与冲激响应的卷积。

$y(t) = x(t) * h(t) = \int_{-\infty}^{\infty} x(\tau) h(t - \tau) \, d\tau$

离散时间情形：

$y[n] = x[n] * h[n] = \sum_{k=-\infty}^{\infty} x[k] h[n - k]$

卷积满足交换律、结合律与分配律。

2. 傅里叶级数 (Fourier Series)

2.1 周期信号的傅里叶级数

周期为 $T_0$ 的信号 $x(t)$ 可以表示为复指数函数之和。基频为 $\omega_0 = 2\pi/T_0$。

$x(t) = \sum_{k=-\infty}^{\infty} c_k e^{jk\omega_0 t}$

傅里叶系数:

$c_k = \frac{1}{T_0} \int_{T_0} x(t) e^{-jk\omega_0 t} \, dt$

2.2 方波(Square Wave)示例

振幅为 A、周期为 T、占空比 50% 的方波的傅里叶系数：

$c_k = \begin{cases} A/2 & k = 0 \\ \frac{A}{k\pi} \sin\left(\frac{k\pi}{2}\right) & k \neq 0 \end{cases}$

只存在奇次谐波，且谐波次数越高振幅越小。吉布斯现象(Gibbs Phenomenon)会出现在不连续点附近。

2.3 帕塞瓦尔定理

$\frac{1}{T_0} \int_{T_0} |x(t)|^2 \, dt = \sum_{k=-\infty}^{\infty} |c_k|^2$

时域的平均功率与频域系数的平方和相等。

3. 傅里叶变换 (Fourier Transform)

3.1 连续傅里叶变换 (CTFT)

非周期信号的傅里叶变换对：

$X(j\omega) = \int_{-\infty}^{\infty} x(t) e^{-j\omega t} \, dt \quad \text{(正变换)}$

$x(t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} X(j\omega) e^{j\omega t} \, d\omega \quad \text{(逆变换)}$

3.2 重要变换对

3.3 傅里叶变换的性质

• 线性: $ax_1(t) + bx_2(t) \leftrightarrow aX_1(j\omega) + bX_2(j\omega)$
• 时移: $x(t-t_0) \leftrightarrow e^{-j\omega t_0} X(j\omega)$
• 频移(调制): $e^{j\omega_0 t} x(t) \leftrightarrow X(j(\omega - \omega_0))$
• 卷积定理: $x(t)*h(t) \leftrightarrow X(j\omega)H(j\omega)$
• 乘法定理: $x(t)h(t) \leftrightarrow \frac{1}{2\pi}X(j\omega)*H(j\omega)$
• 微分: $\frac{dx}{dt} \leftrightarrow j\omega X(j\omega)$
• 帕塞瓦尔: $\int|x(t)|^2 dt = \frac{1}{2\pi}\int|X(j\omega)|^2 d\omega$

3.4 用 Python 计算傅里叶变换

import numpy as np
from scipy.fft import fft, ifft, fftfreq, fftshift
import matplotlib.pyplot as plt

# 生成方波信号
fs = 1000       # 采样频率 (Hz)
T = 1.0         # 信号时长 (秒)
t = np.linspace(0, T, int(fs * T), endpoint=False)

# 方波: 周期 0.1 秒 (10Hz)
from scipy.signal import square
x = square(2 * np.pi * 10 * t)

# 计算 FFT
N = len(x)
X = fft(x)
freqs = fftfreq(N, 1/fs)

# 双边频谱 (归一化)
magnitude = np.abs(fftshift(X)) / N
freq_shifted = fftshift(freqs)

plt.figure(figsize=(12, 4))
plt.subplot(1, 2, 1)
plt.plot(t[:100], x[:100])
plt.title('Square Wave (time domain)')
plt.xlabel('Time (s)')

plt.subplot(1, 2, 2)
plt.plot(freq_shifted, magnitude)
plt.xlim(-100, 100)
plt.title('Fourier Transform (frequency domain)')
plt.xlabel('Frequency (Hz)')
plt.tight_layout()
plt.show()

4. 离散傅里叶变换 (DFT) & FFT

4.1 离散时间傅里叶变换 (DTFT)

离散时间信号 $x[n]$ 的 DTFT:

$X(e^{j\omega}) = \sum_{n=-\infty}^{\infty} x[n] e^{-j\omega n}$

逆变换：

$x[n] = \frac{1}{2\pi} \int_{-\pi}^{\pi} X(e^{j\omega}) e^{j\omega n} \, d\omega$

DTFT 的频率 $\omega$ 是连续的，且以 $2\pi$ 为周期。

4.2 DFT 的定义

$N$ 点 DFT 是把 DTFT 在 $N$ 个频率上做均匀采样得到的结果。

$X[k] = \sum_{n=0}^{N-1} x[n] e^{-j2\pi kn/N}, \quad k = 0, 1, \ldots, N-1$

逆 DFT:

$x[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{j2\pi kn/N}, \quad n = 0, 1, \ldots, N-1$

4.3 FFT 算法

直接计算 DFT 的复杂度是 $O(N^2)$，而 FFT(Cooley-Tukey 算法)是 $O(N \log N)$。当 $N = 2^m$ 时用分治法(Divide and Conquer)实现。

旋转因子(Twiddle Factor): $W_N^k = e^{-j2\pi k/N}$

4.4 频谱分析示例

import numpy as np
from scipy.fft import fft, fftfreq
import matplotlib.pyplot as plt

# 生成信号: 50Hz + 120Hz 复合信号
fs = 1000           # 采样频率 (Hz)
N = fs              # 采样点数 (1秒)
t = np.linspace(0, 1, N, endpoint=False)

# 复合信号 (含噪声)
signal = (np.sin(2 * np.pi * 50 * t)
          + 0.5 * np.sin(2 * np.pi * 120 * t)
          + 0.2 * np.random.randn(N))

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

# 单边频谱 (仅正频率)
half = N // 2
amplitude = 2.0 / N * np.abs(yf[:half])
freq_pos = xf[:half]

plt.figure(figsize=(10, 4))
plt.plot(freq_pos, amplitude)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Amplitude')
plt.title('Single-Sided Amplitude Spectrum')
plt.axvline(x=50, color='r', linestyle='--', label='50 Hz')
plt.axvline(x=120, color='g', linestyle='--', label='120 Hz')
plt.legend()
plt.grid(True)
plt.show()

# 频率分辨率: delta_f = fs / N = 1 Hz
print(f"Frequency resolution: {fs/N} Hz")
print(f"Nyquist frequency: {fs/2} Hz")

4.5 采样定理 (Nyquist-Shannon)

采样频率 $f_s$ 必须至少是信号最大频率 $f_{max}$ 的两倍，才能在没有混叠(Aliasing)的情况下复原信号。

$f_s \geq 2 f_{max} \quad \text{(Nyquist 条件)}$

5. 拉普拉斯变换 (Laplace Transform)

5.1 定义

双边(Bilateral)拉普拉斯变换：

$X(s) = \int_{-\infty}^{\infty} x(t) e^{-st} \, dt, \quad s = \sigma + j\omega$

单边(Unilateral)拉普拉斯变换：

$X(s) = \int_{0^-}^{\infty} x(t) e^{-st} \, dt$

5.2 重要变换对

5.3 传递函数 (Transfer Function)

LTI 系统在 s-域中的输入输出关系：

$H(s) = \frac{Y(s)}{X(s)} = \frac{b_m s^m + b_{m-1} s^{m-1} + \cdots + b_0}{a_n s^n + a_{n-1} s^{n-1} + \cdots + a_0}$

极点(Poles): 使 $H(s) \to \infty$ 的 $s$ 值(分母 = 0)

零点(Zeros): 使 $H(s) = 0$ 的 $s$ 值(分子 = 0)

极点的位置决定系统的稳定性。若所有极点都位于 s-平面左半平面(LHP)，则系统 BIBO 稳定。

5.4 二阶系统标准型

$H(s) = \frac{\omega_n^2}{s^2 + 2\zeta\omega_n s + \omega_n^2}$

• $\omega_n$: 自然频率(Natural Frequency)
• $\zeta$: 阻尼比(Damping Ratio)
  • $\zeta < 1$: 欠阻尼(Underdamped) — 振荡响应
  • $\zeta = 1$: 临界阻尼(Critically Damped) — 最快的无振荡响应
  • $\zeta > 1$: 过阻尼(Overdamped) — 缓慢的无振荡响应

6. Z 变换

6.1 定义

离散时间信号 $x[n]$ 的 Z 变换:

$X(z) = \sum_{n=-\infty}^{\infty} x[n] z^{-n}, \quad z \in \mathbb{C}$

将 $z = e^{j\omega}$ 代入(单位圆上)时，与 DTFT 相同。

6.2 Z 变换与拉普拉斯变换的关系

以采样周期 $T_s$ 对连续时间系统做离散化时：

$z = e^{sT_s}$

由这一关系，s-平面与 z-平面相互对应。

• s-平面的虚轴 $\leftrightarrow$ z-平面的单位圆
• s-平面的左半平面 $\leftrightarrow$ z-平面的单位圆内部

6.3 稳定性分析

离散时间系统要稳定，所有极点都必须位于单位圆内部。

$|z_i| < 1 \quad \text{对所有极点成立}$

6.4 重要 Z 变换对

7. 数字滤波器设计

7.1 FIR vs IIR 滤波器

FIR (Finite Impulse Response) 滤波器

$y[n] = \sum_{k=0}^{M} b_k x[n-k]$

• 冲激响应长度有限
• 始终稳定(无反馈)
• 可实现线性相位
• 要获得陡峭的截止特性，需要较高的阶数

IIR (Infinite Impulse Response) 滤波器

$y[n] = \sum_{k=0}^{M} b_k x[n-k] - \sum_{k=1}^{N} a_k y[n-k]$

• 含反馈，冲激响应无限长
• 用较低阶数即可获得陡峭的截止特性
• 相位非线性(可能产生相位失真)
• 设计不当时可能不稳定

7.2 滤波器类型

• 低通 (LPF): 通过低频，阻断高频
• 高通 (HPF): 通过高频，阻断低频
• 带通 (BPF): 仅通过特定频段
• 带阻 (BEF/Notch): 阻断特定频段

7.3 滤波器设计函数

巴特沃斯(Butterworth): 通带内响应尽可能平坦

切比雪夫(Chebyshev) I 型: 通带内允许波纹，截止更陡峭

椭圆(Elliptic): 通带与阻带都存在波纹，以最低阶数获得最大衰减

7.4 用 scipy.signal 设计滤波器

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

# 采样频率
fs = 1000  # Hz

# 巴特沃斯低通滤波器 (4阶，截止频率 100Hz)
nyq = fs / 2
cutoff = 100 / nyq  # 归一化频率 (0~1)
b_butter, a_butter = signal.butter(N=4, Wn=cutoff, btype='low')

# 切比雪夫 I 型 (4阶，1dB 波纹，截止频率 100Hz)
b_cheby, a_cheby = signal.cheby1(N=4, rp=1, Wn=cutoff, btype='low')

# 计算频率响应
w, h_butter = signal.freqz(b_butter, a_butter, fs=fs)
w, h_cheby = signal.freqz(b_cheby, a_cheby, fs=fs)

plt.figure(figsize=(10, 5))
plt.subplot(1, 2, 1)
plt.plot(w, 20 * np.log10(np.abs(h_butter)), label='Butterworth')
plt.plot(w, 20 * np.log10(np.abs(h_cheby)), label='Chebyshev I')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude (dB)')
plt.title('Filter Frequency Response')
plt.legend()
plt.grid(True)
plt.ylim(-80, 5)

# 对含噪信号进行滤波
t = np.linspace(0, 1, fs, endpoint=False)
clean = np.sin(2 * np.pi * 50 * t)
noisy = clean + 0.5 * np.sin(2 * np.pi * 300 * t)
filtered = signal.filtfilt(b_butter, a_butter, noisy)

plt.subplot(1, 2, 2)
plt.plot(t[:200], noisy[:200], alpha=0.5, label='Noisy')
plt.plot(t[:200], filtered[:200], label='Filtered')
plt.xlabel('Time (s)')
plt.legend()
plt.title('Filtering Result')
plt.tight_layout()
plt.show()

7.5 FIR 滤波器设计 (窗函数法)

from scipy.signal import firwin, freqz

# 使用 Hamming 窗设计的 FIR 低通滤波器
# 抽头数 = 51, 截止频率 = 100Hz
numtaps = 51
cutoff_hz = 100
b_fir = firwin(numtaps, cutoff_hz / (fs/2), window='hamming')

# 频率响应
w_fir, h_fir = freqz(b_fir, [1.0], fs=fs)
print(f"FIR filter order: {numtaps - 1}")
print(f"Linear phase: True (symmetric coefficients)")

8. 控制系统基础

8.1 开环 vs 闭环控制

开环(Open-Loop)控制: 输出不会影响控制动作。结构简单，但对扰动或建模误差较为脆弱。

闭环(Closed-Loop)控制: 测量输出并与参考输入比较，从而使误差最小化。借助反馈(Feedback)获得更高的鲁棒性。

基本反馈回路：

$\frac{Y(s)}{R(s)} = \frac{G(s)}{1 + G(s)H(s)}$

其中 $G(s)$ 是被控对象(Plant)，$H(s)$ 是传感器传递函数。

8.2 稳定性: 劳斯-赫尔维茨判据

若特征方程 $1 + G(s)H(s) = 0$ 的所有根都位于 s-平面左半平面，则系统稳定。

利用劳斯阵列(Routh Array)可仅凭系数判断极点位置。第一列符号变化的次数等于位于右半平面的极点个数。

8.3 波特图 (Bode Plot)

以对数刻度表示频率响应 $H(j\omega)$。

• 增益图: $20\log_{10}|H(j\omega)|$ [dB] vs $\log_{10}\omega$
• 相位图: $\angle H(j\omega)$ [度] vs $\log_{10}\omega$

稳定裕度:

• 相位裕度(Phase Margin, PM): 增益为 0 dB 时，相位距离 -180 度还有多少余量
• 增益裕度(Gain Margin, GM): 相位为 -180 度时，增益距离 0 dB 还有多少余量

通常认为 PM > 45 度、GM > 6 dB 的系统是稳定的。

from scipy import signal
import matplotlib.pyplot as plt
import numpy as np

# 二阶系统传递函数: H(s) = 10 / (s^2 + 2s + 10)
num = [10]
den = [1, 2, 10]
sys = signal.TransferFunction(num, den)

# 波特图
w, mag, phase = signal.bode(sys)

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 6))

ax1.semilogx(w, mag)
ax1.set_ylabel('Magnitude (dB)')
ax1.set_title('Bode Plot')
ax1.grid(True, which='both')

ax2.semilogx(w, phase)
ax2.set_ylabel('Phase (degrees)')
ax2.set_xlabel('Frequency (rad/s)')
ax2.grid(True, which='both')

plt.tight_layout()
plt.show()

9. PID 控制器

9.1 PID 控制动作

PID 控制器是三种动作之和。

$u(t) = K_p e(t) + K_i \int_0^t e(\tau) d\tau + K_d \frac{de(t)}{dt}$

其中 $e(t) = r(t) - y(t)$ 是误差信号。

• 比例(P) 控制: 与当前误差成比例的控制输出。响应快，但会残留稳态误差。
• 积分(I) 控制: 消除累积误差，使稳态误差趋于零。响应会变慢。
• 微分(D) 控制: 对误差变化率作出反应以抑制超调。对噪声敏感。

传递函数:

$C(s) = K_p + \frac{K_i}{s} + K_d s = K_p\left(1 + \frac{1}{T_i s} + T_d s\right)$

9.2 齐格勒-尼科尔斯调参法

方法 1 (开环阶跃响应):

从过程响应曲线求出延迟时间 $L$ 与时间常数 $T$ 后：

方法 2 (闭环极限灵敏度):

令 $K_i = K_d = 0$，逐步增大 $K_p$，求出出现持续振荡的临界增益 $K_u$ 与振荡周期 $T_u$。

9.3 用 Python 模拟 PID

import numpy as np
import matplotlib.pyplot as plt

class PIDController:
    def __init__(self, kp, ki, kd, dt):
        self.kp = kp
        self.ki = ki
        self.kd = kd
        self.dt = dt
        self.integral = 0.0
        self.prev_error = 0.0

    def update(self, setpoint, measured):
        error = setpoint - measured
        self.integral += error * self.dt
        derivative = (error - self.prev_error) / self.dt
        output = self.kp * error + self.ki * self.integral + self.kd * derivative
        self.prev_error = error
        return output

# 一阶被控对象: dy/dt = -y + u (时间常数 1秒)
def plant_step(y, u, dt):
    return y + dt * (-y + u)

# 仿真参数
dt = 0.01
t = np.arange(0, 10, dt)
setpoint = 1.0

# PID 增益设置 (基于齐格勒-尼科尔斯法调参)
pid = PIDController(kp=2.0, ki=1.0, kd=0.5, dt=dt)
y = 0.0
y_list = []

for _ in t:
    u = pid.update(setpoint, y)
    u = np.clip(u, -10, 10)   # 执行器饱和限制
    y = plant_step(y, u, dt)
    y_list.append(y)

plt.figure(figsize=(10, 4))
plt.plot(t, y_list, label='Plant Output')
plt.axhline(setpoint, color='r', linestyle='--', label='Setpoint')
plt.xlabel('Time (s)')
plt.ylabel('Output')
plt.title('PID Control Simulation')
plt.legend()
plt.grid(True)
plt.show()

10. 状态空间表示 (State Space)

10.1 状态方程

连续时间 LTI 系统的状态空间表示:

$\dot{\mathbf{x}}(t) = \mathbf{A}\mathbf{x}(t) + \mathbf{B}\mathbf{u}(t)$ $\mathbf{y}(t) = \mathbf{C}\mathbf{x}(t) + \mathbf{D}\mathbf{u}(t)$

• $\mathbf{x}$: 状态向量 ($n \times 1$)
• $\mathbf{u}$: 输入向量 ($m \times 1$)
• $\mathbf{y}$: 输出向量 ($p \times 1$)
• $\mathbf{A}$: 系统矩阵 ($n \times n$)
• $\mathbf{B}$: 输入矩阵 ($n \times m$)
• $\mathbf{C}$: 输出矩阵 ($p \times n$)
• $\mathbf{D}$: 直接传递矩阵 ($p \times m$)

与传递函数的关系: $H(s) = \mathbf{C}(s\mathbf{I} - \mathbf{A})^{-1}\mathbf{B} + \mathbf{D}$

10.2 可观测性与可控性

可控性(Controllability): 能否在有限时间内，从任意初始状态转移到期望的状态?

若可控性矩阵 $\mathcal{C} = [\mathbf{B} \; \mathbf{AB} \; \mathbf{A}^2\mathbf{B} \; \cdots \; \mathbf{A}^{n-1}\mathbf{B}]$ 的秩为 $n$，则系统完全可控。

可观测性(Observability): 能否仅凭输出确定初始状态?

若可观测性矩阵 $\mathcal{O} = [\mathbf{C}^T \; (\mathbf{CA})^T \; \cdots \; (\mathbf{CA}^{n-1})^T]^T$ 的秩为 $n$，则系统完全可观测。

10.3 线性二次调节器 (LQR)

使代价函数最小化的最优状态反馈控制:

$J = \int_0^{\infty} (\mathbf{x}^T \mathbf{Q} \mathbf{x} + \mathbf{u}^T \mathbf{R} \mathbf{u}) \, dt$

最优控制输入: $\mathbf{u}^* = -\mathbf{K}\mathbf{x}$

$\mathbf{Q}$ 是状态误差权重，$\mathbf{R}$ 是控制输入能量权重。

import numpy as np
from scipy import linalg

# 双积分器系统示例
# x'' = u (二阶系统)
A = np.array([[0, 1], [0, 0]])
B = np.array([[0], [1]])
C = np.array([[1, 0]])
D = np.array([[0]])

# LQR 权重
Q = np.diag([10, 1])  # 更重视位置误差
R = np.array([[1]])

# 求解黎卡提方程
P = linalg.solve_continuous_are(A, B, Q, R)
K = np.linalg.inv(R) @ B.T @ P

print("LQR Gain K:", K)
print("Closed-loop poles:", np.linalg.eigvals(A - B @ K))

11. AI/ML 与信号处理

11.1 用于一维信号的 CNN

一维卷积神经网络非常擅长识别心电图(ECG)、地震信号、语音等时间序列信号中的模式。

import torch
import torch.nn as nn

class Signal1DCNN(nn.Module):
    def __init__(self, num_classes=5):
        super().__init__()
        self.conv_layers = nn.Sequential(
            nn.Conv1d(1, 32, kernel_size=7, padding=3),
            nn.BatchNorm1d(32),
            nn.ReLU(),
            nn.MaxPool1d(2),
            nn.Conv1d(32, 64, kernel_size=5, padding=2),
            nn.BatchNorm1d(64),
            nn.ReLU(),
            nn.MaxPool1d(2),
            nn.Conv1d(64, 128, kernel_size=3, padding=1),
            nn.ReLU(),
            nn.AdaptiveAvgPool1d(8)
        )
        self.classifier = nn.Sequential(
            nn.Flatten(),
            nn.Linear(128 * 8, 256),
            nn.ReLU(),
            nn.Dropout(0.5),
            nn.Linear(256, num_classes)
        )

    def forward(self, x):
        x = self.conv_layers(x)
        return self.classifier(x)

# 使用示例
model = Signal1DCNN(num_classes=5)
# 输入形状: (batch_size, channels=1, sequence_length=1000)
x = torch.randn(8, 1, 1000)
output = model(x)
print("Output shape:", output.shape)  # (8, 5)

11.2 STFT 与 Mel 频谱图

在语音 AI 中，广泛采用将信号转换为时频表示后再交由 CNN 处理的方法。

import librosa
import librosa.display
import matplotlib.pyplot as plt
import numpy as np

# 加载音频
# y, sr = librosa.load('audio.wav', sr=22050)

# STFT (Short-Time Fourier Transform)
# D = librosa.stft(y)
# S_db = librosa.amplitude_to_db(np.abs(D), ref=np.max)

# Mel 频谱图
# mel_spec = librosa.feature.melspectrogram(y=y, sr=sr, n_mels=128)
# mel_db = librosa.power_to_db(mel_spec, ref=np.max)

# MFCC (Mel-Frequency Cepstral Coefficients)
# mfcc = librosa.feature.mfcc(y=y, sr=sr, n_mfcc=13)

print("STFT window: Hann window, hop_length=512, n_fft=2048")
print("Mel-filterbank: 128 mel bands, 0Hz ~ Nyquist")
print("MFCC: 13 coefficients commonly used for speech recognition")

11.3 信号处理 + 深度学习的应用案例

• 心电图(ECG)心律失常分类: 用 1D CNN + LSTM 检测心律异常
• 异常振动检测: 先做 FFT 特征提取，再用自编码器进行异常检测
• 语音识别: MFCC + Transformer (Whisper, wav2vec 2.0)
• 地震检测: 对地震仪信号中的 P 波/S 波自动分类
• 脑电(EEG)分析: 用各频段的功率谱分类专注度/睡眠状态

12. 小测验

参考文献

1. Oppenheim, A. V. & Schafer, R. W. — Discrete-Time Signal Processing (3rd ed.), Prentice Hall, 2010
2. Ogata, K. — Modern Control Engineering (5th ed.), Prentice Hall, 2010
3. Proakis, J. G. & Manolakis, D. G. — Digital Signal Processing, Prentice Hall, 2007
4. SciPy Signal Processing Documentation — https://docs.scipy.org/doc/scipy/reference/signal.html
5. MIT OpenCourseWare 6.003 — Signals and Systems — https://ocw.mit.edu/courses/6-003-signals-and-systems-fall-2011/
6. NumPy FFT 文档 — https://numpy.org/doc/stable/reference/routines.fft.html