别再死记硬背了！用Python+NumPy快速验证三大变换（傅里叶/拉普拉斯/Z变换）公式

用Python+NumPy实战验证信号三大变换：从公式恐惧到可视化理解

记得第一次翻开《信号与系统》教材时，我被傅里叶变换那一页的积分符号和拉普拉斯变换的收敛域条件彻底击垮。直到某天在实验室里，看到学长用三行Python代码生成了完美的频谱图，才恍然大悟——这些看似恐怖的数学工具，本质上都是可以被计算机"计算"和"可视化"的客观现象。本文将带你用NumPy和SciPy，把教材里的公式表变成可交互的代码实验，让变换公式从记忆负担变为直观工具。

1. 环境配置与基础信号生成

工欲善其事，必先利其器。在开始验证变换公式前，我们需要配置好Python科学计算环境。推荐使用Anaconda发行版，它已经集成了我们所需的大部分工具包：

conda create -n signal_analysis python=3.8 numpy scipy matplotlib jupyter
conda activate signal_analysis

接下来在Jupyter Notebook中，我们先定义几个基础信号生成函数：

import numpy as np
import matplotlib.pyplot as plt

def unit_impulse(N, n0=0):
    """生成单位脉冲序列"""
    impulse = np.zeros(N)
    impulse[n0] = 1
    return impulse

def exponential_signal(a, N):
    """生成指数信号a^n"""
    n = np.arange(N)
    return a ** n

def sinusoidal_signal(freq, fs, N):
    """生成正弦信号"""
    t = np.arange(N)/fs
    return np.sin(2*np.pi*freq*t)

注意：采样率fs的选择需满足奈奎斯特准则，应大于信号最高频率的两倍

通过这几个基础函数，我们已经可以构造出教材中80%的示例信号。比如要生成一个衰减指数信号e^(-0.5t)，只需：

t = np.linspace(0, 10, 1000)
signal = np.exp(-0.5 * t)
plt.plot(t, signal)
plt.title('Decaying Exponential Signal')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')

2. 傅里叶变换的数值验证

傅里叶变换将时域信号分解为频率成分，传统手工计算需要处理复杂的积分运算。而使用NumPy的FFT模块，我们可以轻松获得离散傅里叶变换结果。让我们验证几个经典变换对：

案例1：指数衰减信号的傅里叶变换

理论公式：e^(-at) ↔ 1/(a + jω)

a = 0.5
t = np.linspace(0, 10, 1000, endpoint=False)
signal = np.exp(-a * t)

# 计算FFT
fft_result = np.fft.fft(signal)
freqs = np.fft.fftfreq(len(t), t[1]-t[0])

# 理论值
omega = 2 * np.pi * freqs
theoretical = 1 / (a + 1j*omega)

# 可视化对比
plt.figure(figsize=(12,4))
plt.subplot(121)
plt.plot(freqs[:500], np.abs(fft_result)[:500], label='FFT Result')
plt.plot(freqs[:500], np.abs(theoretical)[:500], 'r--', label='Theoretical')
plt.title('Magnitude Spectrum')
plt.legend()

plt.subplot(122)
plt.plot(freqs[:500], np.angle(fft_result)[:500], label='FFT Result')
plt.plot(freqs[:500], np.angle(theoretical)[:500], 'r--', label='Theoretical')
plt.title('Phase Spectrum')
plt.legend()

常见问题排查表 ：

对于更复杂的信号，如矩形窗函数，我们可以观察到著名的sinc函数频谱：

def rect_window(t, T):
    return np.where(np.abs(t) <= T/2, 1, 0)

t = np.linspace(-5, 5, 1000)
signal = rect_window(t, 2)  # 宽度为2的矩形窗

# 计算并绘制频谱
spectrum = np.fft.fftshift(np.fft.fft(np.fft.fftshift(signal)))
freq = np.fft.fftshift(np.fft.fftfreq(len(t), t[1]-t[0]))

plt.plot(freq, np.abs(spectrum))
plt.title('Spectrum of Rectangular Window')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')

3. 拉普拉斯变换的数值逼近

拉普拉斯变换可以视为傅里叶变换的推广，通过引入衰减因子处理不满足绝对可积的信号。虽然NumPy没有直接提供拉普拉斯变换函数，但我们可以通过定义积分核来实现数值计算。

双边指数信号的拉普拉斯变换验证 ：

理论公式：e^(-a|t|) ↔ 2a/(s² - a²)

from scipy.integrate import quad

def laplace_transform(f, s, t_start=0, t_end=np.inf):
    """数值计算单边拉普拉斯变换"""
    integrand = lambda t: f(t) * np.exp(-s * t)
    result, _ = quad(integrand, t_start, t_end)
    return result

# 定义双边指数信号
a = 0.5
f = lambda t: np.exp(-a * np.abs(t))

# 选择s值进行计算
s_values = np.linspace(0.1, 2, 20)
numeric_results = [laplace_transform(f, s) for s in s_values]
theoretical_results = 2*a / (s_values**2 - a**2)

# 绘制对比图
plt.plot(s_values, numeric_results, 'bo', label='Numeric')
plt.plot(s_values, theoretical_results, 'r-', label='Theoretical')
plt.legend()
plt.xlabel('s value')
plt.ylabel('Laplace Transform')
plt.title('Verification of Bilateral Exponential Transform')

对于更实用的场景，我们可以利用拉普拉斯变换求解微分方程的初值问题。例如，考虑一个简单的RC电路：

R = 1e3  # 1kΩ
C = 1e-6  # 1μF
V0 = 5    # 初始电压5V

# 理论解：V(t) = V0 * exp(-t/RC)
t = np.linspace(0, 5e-3, 100)
v_theoretical = V0 * np.exp(-t/(R*C))

# 通过拉普拉斯变换求解
def rc_circuit(t):
    # 数值逆变换近似
    sigma = 1/t[-1]  # 衰减因子
    s = sigma + 1j*np.pi*np.arange(-100,100)/t[-1]
    Fs = V0 / (1 + R*C*s)  # 拉普拉斯域解
    ft = np.fft.ifft(Fs).real * np.exp(sigma*t)
    return ft

v_numeric = rc_circuit(t)

plt.plot(t*1000, v_theoretical, label='Theoretical')
plt.plot(t*1000, v_numeric, 'r--', label='Laplace Solution')
plt.xlabel('Time (ms)')
plt.ylabel('Voltage (V)')
plt.legend()
plt.title('RC Circuit Response via Laplace Transform')

4. Z变换与离散系统分析

Z变换是分析离散时间系统的强大工具。让我们用Python验证几个基本序列的Z变换关系。

案例：指数序列的Z变换

理论公式：a^n u[n] ↔ z/(z - a)

def z_transform(sequence, z):
    """计算序列在给定z点的Z变换"""
    n = np.arange(len(sequence))
    return np.sum(sequence * (z ** (-n)))

a = 0.8
N = 50
n = np.arange(N)
x = a ** n  # 指数序列

# 选择z值进行评估
z_magnitude = 1.0  # 单位圆上
z_angles = np.linspace(0, 2*np.pi, 200)
z_values = z_magnitude * np.exp(1j * z_angles)

# 计算数值Z变换和理论值
X_numeric = np.array([z_transform(x, z) for z in z_values])
X_theoretical = z_values / (z_values - a)

plt.figure(figsize=(12,5))
plt.subplot(121)
plt.plot(z_angles, np.abs(X_numeric), label='Numeric')
plt.plot(z_angles, np.abs(X_theoretical), 'r--', label='Theoretical')
plt.title('Magnitude Response')
plt.xlabel('Frequency (rad)')
plt.legend()

plt.subplot(122)
plt.plot(z_angles, np.angle(X_numeric), label='Numeric')
plt.plot(z_angles, np.angle(X_theoretical), 'r--', label='Theoretical')
plt.title('Phase Response')
plt.xlabel('Frequency (rad)')
plt.legend()

Z变换应用：系统稳定性分析

通过计算系统函数的极点位置，可以判断系统稳定性：

def plot_pole_zero(b, a):
    """绘制系统函数的零极点图"""
    zeros = np.roots(b)
    poles = np.roots(a)
    
    plt.figure(figsize=(6,6))
    unit_circle = plt.Circle((0,0), 1, fill=False, linestyle='--')
    plt.gca().add_patch(unit_circle)
    
    plt.scatter(np.real(zeros), np.imag(zeros), marker='o', facecolors='none', edgecolors='b', label='Zeros')
    plt.scatter(np.real(poles), np.imag(poles), marker='x', color='r', label='Poles')
    plt.legend()
    plt.axis('equal')
    plt.grid()
    plt.title('Pole-Zero Plot')

# 示例系统：y[n] = 0.9y[n-1] + x[n]
b = [1]    # 分子系数
a = [1, -0.9]  # 分母系数
plot_pole_zero(b, a)

5. 综合应用：音频信号频谱分析

将三大变换知识综合应用，我们可以构建一个简单的音频频谱分析工具。以下代码演示如何分析录音信号的时频特性：

import sounddevice as sd
from scipy.io import wavfile

# 录制或加载音频文件
fs, audio = wavfile.read('sample.wav')  # 或使用sd.rec()录制
audio = audio[:,0] if audio.ndim > 1 else audio  # 取单声道

# 短时傅里叶变换分析
window_size = 1024
hop_length = 512
n_windows = (len(audio) - window_size) // hop_length + 1

spectrogram = np.zeros((window_size//2 + 1, n_windows))
for i in range(n_windows):
    segment = audio[i*hop_length : i*hop_length+window_size]
    windowed = segment * np.hamming(window_size)
    spectrum = np.abs(np.fft.rfft(windowed))
    spectrogram[:,i] = 20 * np.log10(spectrum + 1e-12)

# 绘制声谱图
plt.figure(figsize=(12,6))
plt.imshow(spectrogram, aspect='auto', origin='lower',
           extent=[0, len(audio)/fs, 0, fs/2], cmap='inferno')
plt.colorbar(label='Intensity (dB)')
plt.xlabel('Time (s)')
plt.ylabel('Frequency (Hz)')
plt.title('Audio Spectrogram')

音频处理常用变换对比表 ：

在完成这些实验后，我习惯将常用的验证代码封装成可复用的函数模块。比如创建一个signal_utils.py文件，包含各种变换的验证函数，这样下次遇到新的公式时，只需几行调用就能快速验证其正确性。这种"实验式学习"方法彻底改变了我对信号处理的理解方式——公式不再是需要死记硬背的符号，而是可以通过代码"触摸"和"观察"的活工具。