傅里叶变换、FFT、S变换与Z变换可视化

本教程通过Python代码实现核心公式+动态图演示，帮助你直观理解这些变换的本质。

一、傅里叶变换（Fourier Transform）

核心思想

将信号分解为不同频率的正弦波分量，揭示信号的频率组成。

连续傅里叶变换公式：

$F(\omega) = \int_{-\infty}^{\infty} f(t)e^{-j\omega t}dt$

Python实现：

python
展开代码
import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, fftfreq
from pylab import mpl
mpl.rcParams['font.sans-serif'] = ['FangSong']  # 指定默认字体
mpl.rcParams['axes.unicode_minus'] = False  # 解决保存图像是负号'-'显示为方块的问题

# 生成合成信号
fs = 1000            # 采样率
T = 1/fs             # 采样周期
L = 1                # 信号长度（秒）
t = np.arange(0, L, T)

# 创建包含2个频率成分的信号
f1, f2 = 50, 120     # 50Hz和120Hz
y = 0.7*np.sin(2*np.pi*f1*t) + np.sin(2*np.pi*f2*t)

# 执行FFT
N = len(y)
yf = fft(y)
xf = fftfreq(N, T)[:N//2]

# 绘制时域和频域图
plt.figure(figsize=(12,6))
plt.subplot(2,1,1)
plt.plot(t, y)
plt.title('时域信号')
plt.xlabel('时间 [s]')

plt.subplot(2,1,2)
plt.plot(xf, 2.0/N * np.abs(yf[0:N//2]))
plt.title('频域信号（FFT）')
plt.xlabel('频率 [Hz]')
plt.tight_layout()
plt.show()

关键观察：

• 频谱图中50Hz和120Hz处出现明显峰值
• 幅值大小反映不同频率成分的能量强度

二、快速傅里叶变换（FFT）

核心优势：

将DFT的复杂度从O(N²)优化到O(N logN)，使实时信号处理成为可能

FFT分辨率实验：

python
展开代码
import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, fftfreq, rfft, rfftfreq
from pylab import mpl

mpl.rcParams['font.sans-serif'] = ['FangSong']  # 指定默认字体
mpl.rcParams['axes.unicode_minus'] = False  # 解决保存图像是负号'-'显示为方块的问题'


def generate_signal(fs=1000, L=1):
    """生成合成信号"""
    T = 1 / fs  # 采样周期
    t = np.arange(0, L, T)  # 时间向量

    # 创建包含2个频率成分的信号 (50Hz 和 120Hz)
    f1, f2 = 50, 120
    y = 0.7 * np.sin(2 * np.pi * f1 * t) + np.sin(2 * np.pi * f2 * t)

    # 添加一些随机噪声
    y += np.random.normal(0, 0.3, len(t))

    return t, y


def perform_fft(t, y):
    """执行FFT分析"""
    N = len(y)
    T = t[1] - t[0]

    # 计算FFT
    yf = fft(y)

    # 获取对应的频率
    xf = fftfreq(N, T)[:N // 2]

    # 计算单边幅值
    amp = 2.0 / N * np.abs(yf[:N // 2])

    return xf, amp


def perform_real_fft(t, y):
    """执行实数FFT优化计算"""
    N = len(y)
    T = t[1] - t[0]

    # 计算实数FFT (针对实数输入的优化版本)
    yf = rfft(y)

    # 获取对应的频率
    xf = rfftfreq(N, T)

    # 计算幅值
    amp = 2.0 / N * np.abs(yf)

    return xf, amp


def plot_results(t, y, xf, amp, title='FFT 分析'):
    """绘制时域和频域图形"""
    plt.figure(figsize=(12, 8))

    # 时域信号
    plt.subplot(2, 1, 1)
    plt.plot(t, y)
    plt.title('时域信号')
    plt.xlabel('时间 [s]')
    plt.ylabel('幅值')
    plt.grid(True)

    # 频域信号
    plt.subplot(2, 1, 2)
    plt.plot(xf, amp)
    plt.title(title)
    plt.xlabel('频率 [Hz]')
    plt.ylabel('幅值')
    plt.grid(True)
    plt.tight_layout()
    plt.show()


def calculate_psd(xf, amp):
    """计算功率谱密度"""
    power = amp ** 2 / len(amp)
    return xf, power


def zero_padding_analysis(t, y, pad_factor=2):
    """零填充分析"""
    # 零填充信号
    y_padded = np.pad(y, (0, int(len(y) * (pad_factor - 1))), 'constant')  # 修正这里
    t_padded = np.linspace(t[0], t[-1], len(y_padded))  # 修正这里

    # 计算FFT
    N = len(y_padded)
    T = t_padded[1] - t_padded[0]
    yf_padded = rfft(y_padded)
    xf_padded = rfftfreq(N, T)
    amp_padded = 2.0 / N * np.abs(yf_padded)

    return t_padded, xf_padded, amp_padded, y_padded  # 返回y_padded用于绘图


def main():
    # 参数设置
    fs = 1000  # 采样率
    L = 1  # 信号长度（秒）

    # 生成信号
    t, y = generate_signal(fs, L)

    # 执行FFT
    xf, amp = perform_real_fft(t, y)

    # 绘制基本分析结果
    plot_results(t, y, xf, amp, 'FFT 分析 (原始信号)')

    # 显示主要频率成分
    print("\n检测到的主要频率成分：")
    peaks = np.where(amp > 0.3)[0]  # 幅值大于0.3的峰
    for i in peaks:
        print(f"频率: {xf[i]:.2f} Hz, 幅值: {amp[i]:.4f}")

    # 功率谱密度分析
    freq_psd, power = calculate_psd(xf, amp)
    plot_results(t, y, freq_psd, power, '功率谱密度')

    # 零填充分析
    t_padded, xf_padded, amp_padded, y_padded = zero_padding_analysis(t, y, pad_factor=5)  # 接收y_padded
    plot_results(t_padded, y_padded, xf_padded, amp_padded, '零填充后的FFT分析')

    # 对比图
    plt.figure(figsize=(12, 6))
    plt.plot(xf, amp, label='原始FFT')
    plt.plot(xf_padded, amp_padded, label='零填充FFT')
    plt.title('原始信号与零填充信号FFT对比')
    plt.xlabel('频率 [Hz]')
    plt.ylabel('幅值')
    plt.legend()
    plt.grid()
    plt.show()


if __name__ == "__main__":
    main()

实验结论：

• 零填充（Zero-padding）可提高频谱显示分辨率
• 但不会增加真实频率分辨率（受原始信号长度限制）

展开代码
检测到的主要频率成分：
频率: 50.00 Hz, 幅值: 0.6987
频率: 120.00 Hz, 幅值: 1.0120

三、S变换（Stockwell Transform）

时频分析利器：

结合短时傅里叶变换（STFT）和小波变换的优点，具有可逆性

Python实现：

python
展开代码
from scipy.signal import stft

# 生成非平稳信号（频率随时间变化）
t = np.linspace(0, 10, 1000)
f = np.linspace(1, 50, len(t))
y = np.sin(2*np.pi*f*t)

# 计算S变换
frequencies, times, Zxx = stft(y, fs=1/(t[1]-t[0]), nperseg=128)

# 绘制时频图
plt.figure(figsize=(10,6))
plt.pcolormesh(times, frequencies, np.abs(Zxx), shading='gouraud')
plt.title('S变换时频图')
plt.ylabel('频率 [Hz]')
plt.xlabel('时间 [sec]')
plt.colorbar(label='幅值')
plt.show()

关键特性：

• 横轴：时间变化
• 纵轴：频率分布
• 颜色：能量强度
• 能清晰捕捉频率随时间的变化趋势

四、Z变换（Z-Transform）

离散系统分析工具：

将差分方程转换为代数方程，分析数字滤波器的稳定性

零极点图可视化：

python
展开代码
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import tf2zpk

# 定义系统函数系数
b = [1, -1.5]    # 分子系数（零点）
a = [1, -0.9, 0.8]  # 分母系数（极点）

# 计算零极点
z, p, k = tf2zpk(b, a)

# 绘制零极点图
plt.figure(figsize=(6,6))
plt.plot(np.real(z), np.imag(z), 'o', markersize=10, markerfacecolor='none', markeredgecolor='r')
plt.plot(np.real(p), np.imag(p), 'x', markersize=10, markeredgecolor='b')

# 绘制单位圆
unit_circle = plt.Circle((0,0), 1, color='black', fill=False)
plt.gca().add_patch(unit_circle)
plt.axhline(0, color='black')
plt.axvline(0, color='black')
plt.grid(True)
plt.title('零极点图')
plt.xlabel('实部')
plt.ylabel('虚部')
plt.axis('equal')
plt.show()

稳定性判断：

• 极点（×）全部在单位圆内 → 系统稳定
• 零点（○）位置影响频率响应特性

五、变换方法对比总结

变换类型	适用领域	时频分辨率	稳定性分析	可逆性
傅里叶变换	平稳信号分析	无时间分辨率	否	是
FFT	快速频谱计算	固定分辨率	否	是
S变换	非平稳信号分析	自适应分辨率	否	是
Z变换	离散系统设计	无	是	是

六、扩展练习建议

1. 尝试改变信号的采样率（fs），观察对频谱的影响
2. 添加不同类型的噪声（高斯白噪声、工频干扰等）
3. 设计简单FIR滤波器在频域进行滤波操作
4. 使用S变换分析真实心电信号的时频特征
5. 通过Z变换设计IIR数字滤波器

提示：可使用scipy.signal库中的find_peaks检测频谱峰值，spectrogram生成更丰富的时频图