3.3 FFT 频谱观测实战

3.3 FFT 频谱观测实战

本节摘要：本节给出从 I/Q 样本到频谱图与瀑布图的完整可复用脚本：复数 FFT 加 fftshift 的标准流程、平均化压噪、窗函数选择、瀑布图绘制，并用 FM 广播频段与 433 MHz 遥控信号两个实例演示"频谱读图"的方法。这是全书被引用次数最多的一节。

频谱脚本的最小骨架

第 2 章说过复数样本要成对组合，第 3.1 节讲了频率刻度的算法，现在把它们拼成标准流程。下面这份脚本值得原样保存，后续章节的实验都从它出发：

import numpy as np
import matplotlib.pyplot as plt

def load_iq(fname):
    raw = np.fromfile(fname, dtype=np.uint8).astype(np.float32) - 127.5
    return (raw[0::2] + 1j * raw[1::2]) / 128.0

def psd(iq, fs, nfft=4096, avg=16):
    step = nfft * avg // 4                  # 重叠分段，兼顾平滑与速度
    frames = [iq[i:i+nfft] for i in range(0, len(iq)-nfft, step)][:avg]
    win = np.hanning(nfft)
    spec = np.zeros(nfft)
    for fr in frames:
        spec += np.abs(np.fft.fftshift(np.fft.fft(fr * win))) ** 2
    spec /= len(frames) * np.sum(win ** 2) / fs * 4 / fs
    f = np.fft.fftshift(np.fft.fftfreq(nfft, 1 / fs))
    return f, 10 * np.log10(spec + 1e-15)

iq = load_iq("fm_band.iq")                   # 第 2 章采集的 98 MHz 附近样本
f, p = psd(iq, 2.048e6)
plt.plot(f / 1e3, p)
plt.xlabel("相对中心频率 kHz"); plt.ylabel("相对功率 dB")
plt.title("调频广播频段频谱"); plt.grid(True); plt.show()
​

三个设计点说明：fftshift 把零频移到画面中央（复数频谱天然从负频率开始，不 shift 看图会很别扭）；分段平均是压低噪声随机起伏的关键，单帧频谱抖动剧烈，平均十几帧后噪声底就平得能看清细节；窗函数乘在 FFT 之前，下一小节展开。横轴用"相对中心频率"是 SDR 频谱图的惯例——画面中心永远是你的调谐频率。

窗函数：频谱泄漏与它的药方

FFT 默认假设你的片段是无限周期信号的一小段，但真实信号不会在片段边界恰好归零，边界的不连续会在频谱上拖出一条"裙边"，把强信号的能量洒到邻近频点上——这就是频谱泄漏。不加窗时，一个强台的泄漏裙边足以盖过旁边的弱信号。加窗就是给片段两端乘上渐变到零的权重，消除边界不连续；代价是主瓣变宽、频率分辨率略降。

常用窗的取舍：矩形窗（不加窗）分辨率最高但泄漏最重，仅适合间隔很开的强信号；汉宁窗（hanning）泄漏压制好、主瓣居中，是频谱观测的默认选择；汉明窗（hamming）旁瓣更低一档；布莱克曼窗（blackman）泄漏最小但主瓣最宽，适合"强信号旁边找弱信号"的场景。经验法则：一般观测用汉宁，找深埋信号换布莱克曼。

图：矩形窗与汉宁窗的频谱泄漏对比示意

实例一：FM 广播频段的读图

对 98 MHz 附近采集的样本跑上面的脚本，典型画面是：一条平坦的噪声底上矗立着两三座"平顶山"，每座宽约 180 到 200 kHz——这就是宽带调频信号的频域指纹，平顶来自调频的恒包络特性，宽度由最大频偏决定（75 kHz 频偏的广播信号，卡森带宽约 180 kHz）。读图要点有三：峰的宽度告诉你调制类型；峰顶是否平直告诉你是否过调或过增益；中央那根细尖的直流尖峰是零中频架构的本振泄漏，观测时忽略它或轻微偏移中心频率绕开它。

实例二：433 MHz 遥控信号的瀑布图

瀑布图把频谱沿时间堆叠成图像，横轴频率、纵轴时间、亮度为强度，特别适合捕捉间歇信号。最常见的入门目标是无钥匙遥控钥匙或 433 MHz 温湿度传感器：平时频谱一片死寂，按下按键的瞬间一道亮带闪现。

def waterfall(iq, fs, nfft=1024, rows=200):
    step = (len(iq) - nfft) // rows
    img = np.zeros((rows, nfft))
    win = np.hanning(nfft)
    for r in range(rows):
        seg = iq[r * step : r * step + nfft] * win
        img[r] = 10 * np.log10(np.abs(np.fft.fftshift(np.fft.fft(seg))) ** 2 + 1e-12)
    plt.imshow(img, aspect="auto", cmap="viridis",
               extent=[-fs/2/1e3, fs/2/1e3, rows, 0])
    plt.xlabel("相对频率 kHz"); plt.ylabel("帧序号")
    plt.title("瀑布图"); plt.show()

iq433 = load_iq("remote_433.iq")     # 对 433.92 MHz 采集，期间按几次遥控器
waterfall(iq433, 2.048e6)
​

画面解读：遥控信号通常是 ASK 或 FSK，表现为一条短促的亮带（ASK 单峰）或一按一亮的双线（FSK 两个频点交替），持续时间几十毫秒、间隔数百毫秒重复。从瀑布图里量出"亮带的宽度与持续时间"，就得到了下一步解码所需的全部参数——信号带宽定滤波器宽度，脉冲节奏定解码时钟。第 8 章的卫星 APT 案例会再次用到同一套观测方法。

频谱观测的常见误区

三个高频误区值得点名。其一，把直流尖峰当成强信号：它是零中频的本振泄漏，永远在正中央，换个中心频率它还跟着你。其二，把过增益的平顶削波当成信号很稳：削波在频谱上表现为突然冒出的谐音符，降低增益即消失。其三，把 FFT 点数当分辨率调大就完事：点数提高分辨率的同时，瀑布图时间轴也变慢，捕捉短促信号反而要用小点数。三误区的共同解法是回到定义：频谱图的每个像素都是"某频率、某时刻、某带宽内"的功率，搞清三个坐标的含义，图就不会读错。

常见问题

问题：为什么我画的频谱噪声底一直在跳？

单帧 FFT 的噪声估计方差极大，每帧之间起伏 10 dB 以上都属正常。解决方法就是脚本里的分段平均：十几帧平均后噪声底趋于平稳。若要更稳，把 avg 参数加大到 32 或 64，代价是画面刷新变慢。

问题：fftshift 到底做了什么？

复数信号的 FFT 输出排列是"正频率在前、负频率在后"，直接画图会出现信号被切成两半的画面。fftshift 把零频挪到数组中央，负频率在左、正频率在右，符合人眼阅读习惯。它只重排显示顺序，不改任何数值。

问题：观测到的信号强度单位是什么？

本节脚本的纵轴是相对功率（未做绝对校准），数值只用于比较相对强弱。想要绝对的 dBm 需要用已知功率的标准信号源标定——对多数实验并无必要，相对比较已经够用：比较两个台的强弱、观察增益改变前后的变化，相对值都能给出正确结论。

本节要点回顾

• 标准流程：成对组合 I/Q，加窗，FFT，fftshift，多帧平均，横轴用相对频率。
• 窗函数取舍：汉宁是默认，布莱克曼对付强邻道，矩形窗只适合间隔大的强信号。
• FM 广播指纹：约 180 到 200 kHz 的平顶峰，卡森带宽由 75 kHz 频偏决定。
• 瀑布图看间歇信号：亮带的宽度与持续时间就是解码参数的来源。

频谱告诉你信号"在哪、多宽、什么形状"，下一节讲信息是怎么被"装进"这些形状的——调制的家族图谱。