用Matlab进行无线电信号逆向实战2——立体声 FM 广播的分离与解密 从频谱迷宫到相干解调的避坑指南

在上一期的 Q01 中，我们成功解调了单声道 FM 广播。但现实世界中的无线电信号往往更为复杂，为了在一个载波上传输更多信息，频分复用（FDM）技术应运而生。本期（Q02）我们将踏入立体声 FM 广播的领域。
本题在 SDR CTF 体系中属于模拟广播层的进阶，难度 ★2。它不仅引入了频分复用和相干解调的概念，更是一个绝佳的“信号诊断”实战案例——当实际跑出的频谱与教科书理论严重不符时，我们该如何由表及里地排查问题，而不是对着代码怀疑人生。

一、 题目设定

在 CTF 竞赛或真实的无线电逆向中，拿到一个未知的 IQ 文件，我们不应主观臆断它的类型。以下是本题的题干设定：

• 截获文件：Q02_stereo_fm.iq（复基带信号，叠加了 SNR=25dB 的 AWGN 噪声）下载地址。
• 已知参数：初始告知采样率Fs = 240 kHz（有个小陷阱），信号为恒包络调制，载波已下变频至基带。
• 挑战目标：成功分离出立体声的左声道(L)和右声道®。
• Flag 载体：Flag 隐藏在右声道®音频中，以 800 Hz 载波承载的简化摩斯电码（ON/OFF 序列[1 0 1 1 1 0 1 1 1 0 0 0]）形式发送。
• 预期 Flag：FLAG{FM_STEREO_R}

二、 盲解第一步：频谱“指纹”失效与诊断

面对未知的 IQ 文件，标准的盲分析流程是“时域看包络 → 频域看轮廓 → 鉴频转基带 → 功率谱看指纹”。

1. 常规操作与意外的频谱

读取文件后，我们先观察时域包络，发现幅度恒定，确认是恒包络调制（FM/FSK/PSK）。

于是顺理成章地进行 FM 鉴频（共轭相乘法），提取基带信号，并使用 Welch 法绘制平滑功率谱：

% 读取 IQ 文件fid=fopen('Q02_stereo_fm.iq','r');raw=fread(fid,'float32');fclose(fid);iq=complex(raw(1:2:end),raw(2:2:end));Fs=240e3;% 初始给定的采样率% FM 鉴频 (共轭相乘法)iq=iq-mean(iq);% 去除直流偏置fm_demod=angle(iq(2:end).*conj(iq(1:end-1)));baseband=fm_demod/(2*pi)*Fs;% 转换为 Hz 量级baseband=baseband-mean(baseband);% 去除鉴频器输出的残余直流% 绘制功率谱密度 (Welch 法)[pxx,f]=pwelch(baseband,hann(4096),2048,4096,Fs);figure;plot(f/1e3,10*log10(pxx));xlabel('kHz');ylabel('dB/Hz');title('Fs=240kHz 时的基带功率谱');

按照理论，鉴频后的立体声基带频谱应该呈现标准的三段式“指纹”：0~15kHz 的音频、19kHz 的孤立尖峰、23~53kHz 的对称能量块。

2. 奇怪图像分析：中高频段爆发与采样率混叠

图像表现：

频谱图呈现出“两头高，中间低”的怪异特征：

1. 极低频区域 (0 - 5 kHz)：出现极高能量的尖峰群，瞬间达到约 60 dB/Hz。
2. 中频区域 (5 - 65 kHz)：看似平坦，能量极低，完全没有导频和副载波的影子。
3. 中高频爆发 (70 - 85 kHz)：在76 kHz 附近出现了极强的高能“爆发”，峰值超过 60 dB/Hz，几乎与低频尖峰持平，随后迅速衰减。
   由表及里的根因剖析：
   为什么理论与实测对不上？核心问题出在采样率不足导致的频谱混叠。

• 理论带宽：FM 信号的带宽遵循卡森公式B = 2 ( Δ f + f m ) B = 2(\Delta f + f_m)B=2(Δf+fm​)。本题最大频偏Δ f = 75 \Delta f = 75Δf=75kHz，复合基带最高频率f m = 53 f_m = 53fm​=53kHz。实际带宽约2 × ( 75 + 53 ) = 256 2 \times (75 + 53) = 2562×(75+53)=256kHz。
• 混叠发生：初始采样率 240 kHz低于 256 kHz 的信号带宽。在数字域生成信号时，超出± 120 \pm 120±120kHz（奈奎斯特频率）的频谱分量被折叠回来。
• 非线性恶化：混叠破坏了恒包络特性。当混叠信号进入鉴频器（非线性操作）时，原本折叠的高频分量与 38kHz 副载波发生交调，在 76kHz 产生了极强的杂散峰，彻底掩盖了微弱的 19kHz 导频。
  解决方案：必须提升采样率。我们将发送端与接收端的采样率统一下调至原设计的合理值（例如重新提供F s = 500 Fs = 500Fs=500kHz 的数据，或按比例重采样），确保奈奎斯特带宽覆盖所有信号分量。
  真实数据下载地址。

三、 盲解第二步：修正采样率与发现真实频谱

1. 修正代码与图像表现

在 500 kHz 采样率下，奈奎斯特带宽达到 250 kHz，彻底消除了混叠。再次运行上述pwelch代码，频谱基底噪声变得平稳（约 35 dB/Hz），真实的信号结构显现出来：

1. 极低频 (0-5 kHz) 强信号 (54-56 dB/Hz)：这是完全正常的 (L+R) 主信道，因为左右声道的音频（440Hz, 880Hz, 800Hz）都集中于此。
2. 核心峰值 (36.5 kHz 附近, 55 dB/Hz)：这是 (L-R) 副信道（38 kHz DSB-SC 调制）的能量。由于左右声道差异大且 FM 调制的非线性交调，能量在此处高度集中。
3. 15 - 20 kHz 的轻微隆起 (43 dB/Hz)：这是 19 kHz 导频。

2. 奇怪图像分析：导频为何变成了“小山包”？

图像表现：
为什么 19 kHz 导频在频谱上只是一个微弱的“隆起”，而不是教科书里孤立高耸的“尖峰”？如果仅凭这张图，很容易误判该信号不是立体声。
由表及里的根因剖析：

• 三角噪声效应：FM 鉴频器有一个固有的物理特性——输出噪声功率谱密度与频率的平方成正比（即高频噪声远大于低频噪声）。
• 导频幅度极小：在发送端，导频的幅度被特意设为复合基带最大幅度的 10% 左右，以节省频偏资源。
  在 19 kHz 这个相对较高的频率上，噪声底已经被三角噪声抬高，微弱的导频信号在 Welch 功率谱上就被“糊”成了一个隆起的包，而不是一根清晰的针。
  解决方案：既然直接看频谱看不清，我们就不能仅凭肉眼看频谱，必须用极窄的带通滤波器去把导频“挖”出来作为证据。

四、 盲解第三步：用带通滤波验证隐藏的导频

既然 19kHz 导频在频谱上不明显，我们设计一个极窄的零相位带通滤波器，将其从噪声中强制提取出来，并观察时域波形。

% 设计 19kHz 窄带带通滤波器 (零相位)bpFilt_19k=designfilt('bandpassfir','FilterOrder',500,...'CutoffFrequency1',18.8e3,'CutoffFrequency2',19.2e3,'SampleRate',Fs);pilot_19k=filtfilt(bpFilt_19k,baseband);% 必须用 filtfilt 保证零相位% 观察导频时域波形figure;t_ms=(0:4999)/Fs*1000;% 取前 10ms 数据plot(t_ms,pilot_19k(1:5000));xlabel('时间');ylabel('幅度');title('提取出的 19kHz 导频时域波形');

奇怪图像分析：时域波形的“视觉陷阱”

图像表现：
运行代码后，你会看到一张布满密集振荡线条的图，幅度在 -6000 到 6000 之间剧烈波动，看起来像是一个“实心的蓝色填充区域”。波形在某些时间点达到最大振幅（如 X≈9.8 处峰值 5800），在另一些时间点振幅变得很小甚至接近零（如 X≈7 处幅度收缩）。这看起来像是导频受到了低频调制或发生了拍频！
由表及里的根因剖析：
这是一种典型的视觉错觉。19kHz 是高频信号，每毫秒包含 19 个完整周期。当我们在一个较宽的时间窗口（如 10ms）内绘制几百个周期时，波形密度极高，相邻周期的波峰和波谷在屏幕像素级渲染时会相互重叠，产生类似摩尔纹的干涉现象。

• 图中看到的“包络上下边界”，其实就是正弦波本身的等幅上下限在像素上的压缩表现。
• 看似“剧烈波动的包络”，实际上是高频正弦波在渲染时的混叠视觉效应。19kHz 导频在发送端是未调制的单频正弦波，理论上是严格等幅的。
  解决方案：缩短时间窗口验证。
  若想验证其稳定性，只需将时间轴缩短到 1ms 以内（约 19 个周期），就能看到完美的等幅正弦波。

figure;t_ms_short=(0:499)/Fs*1000;% 取前 1ms 数据plot(t_ms_short,pilot_19k(1:500));xlabel('时间');ylabel('幅度');title('19kHz 导频时域波形（短时间窗口验证）');

验证成功！在短窗口下（代码中窗口还是有些大了），导频呈现出完美的等幅正弦波。至此，“标准 FM 立体声广播”的判定才有了充分的局部证据。

五、 解调实现中的三大“致命陷阱”

判明信号类型后，进入解调阶段。这里同样有三个易错点需要严格规避，否则前功尽弃。

陷阱一：数学相位的暗度陈仓（sin 与 cos 之争）

现象：接收端提取 19 kHz 导频，平方后得到 38 kHz 副载波，将其与 DSB-SC 信号相乘并低通滤波，结果 (L-R) 信号全是 0。
剖析：如果发射端导频和副载波使用了sin：

• 提取导频并平方：P² = sin²(2π·19k·t) = 0.5 - 0.5·cos(2π·38k·t)
• 38 kHz 带通滤波后，本地副载波变成了：-cos(2π·38k·t)
• 发射端的 DSB-SC 信号是：(L-R) · sin(2π·38k·t)
• 相乘解调：[(L-R) · sin(38k)] · [-cos(38k)] = -0.5 · (L-R) · sin(76k)
  乘积变为高频项（76kHz），被低通滤波器彻底滤除！
  避坑指南：在涉及倍频和相干解调的系统中，导频和副载波必须全部使用cos。平方后得到同相cos(38k)，相乘后产生直流分量0.5(L-R)，才能被低通滤波器保留。

陷阱二：滤波器的群延迟失配

现象：修正了 sin/cos 问题后，解调出来的 L-R 依然是一片噪声，摩斯电码能量检测全为 0。
剖析：FIR 滤波器具有群延迟，阶数为N NN的滤波器会使信号延迟N / 2 N/2N/2个采样点。

• DSB-SC 信号经过 23-53kHz 带通滤波，延迟了N 1 / 2 N_1/2N1​/2。
• 导频经过 19kHz 带通滤波，平方后再次经过 38kHz 带通滤波，累计延迟了N 2 / 2 + N 3 / 2 N_2/2 + N_3/2N2​/2+N3​/2。
  两路信号在时间轴上错位。高频周期信号即使错开几个采样点，也会导致巨大的相位差（如从0 ∘ 0^\circ0∘变成90 ∘ 90^\circ90∘），相干解调输出幅度趋近于 0。
  避坑指南：绝不能使用单向的filter函数，必须使用filtfilt（零相位滤波）。正向和反向各滤波一次，相互抵消群延迟，保证相位严格对齐。

陷阱三：硬件声卡的“降维打击”（采样率不支持）

现象：解调成功后，用sound(R, Fs)播放，报错Device Error: Unanticipated host error。
剖析：500 kHz 的采样率远超普通 PC 声卡支持范围（通常最高 48kHz 或 192kHz），底层音频 API 拒绝接受过高采样率的音频流。
避坑指南：输出到声卡前必须降采样。使用resample(R, 48000, Fs)将其转换为 48 kHz，即可顺利播放。

六、 核心解调逻辑演示

基于上述避坑经验，我们提取 L-R 信号并进行立体声解码的核心逻辑如下（非完整脚本，仅展示关键环节）：

% 1. 设计零相位滤波器 (此处省略 designfilt 细节)% lpFilt: 15kHz 低通 | bpFilt_19k | bpFilt_38k | bpFilt_23_53% 2. 提取主信道 (L+R)L_plus_R=filtfilt(lpFilt,baseband);% 3. 提取 19kHz 导频并倍频恢复 38kHz 副载波pilot_19k=filtfilt(bpFilt_19k,baseband);subcarrier_38k=pilot_19k.*pilot_19k;subcarrier_38k=filtfilt(bpFilt_38k,subcarrier_38k)*2;% 乘2补偿幅度% 4. 相干解调副信道dsb_sc=filtfilt(bpFilt_23_53,baseband);L_minus_R=filtfilt(lpFilt,dsb_sc.*subcarrier_38k);% 5. 立体声矩阵解码L=(L_plus_R+L_minus_R)/2;R=(L_plus_R-L_minus_R)/2;% 6. 针对 R 通道进行 800Hz 摩斯电码能量检测% 对 R 进行 800Hz 带通滤波，按符号长度计算短时能量并阈值判决% 最终检测到序列 [1 0 1 1 1 0 1 1 1 0 0 0]% 解出 Flag: FLAG{FM_STEREO_R}

完整代码如下：

% Q02_stereo_fm_rx.m - 解题脚本（零相位滤波 + 音频降采样）clear;clc;% 1. 读取 IQ 文件fid=fopen('Q02_stereo_fm.iq','r');raw=fread(fid,'float32');fclose(fid);iq_data=complex(raw(1:2:end),raw(2:2:end));Fs=500e3;% 采样率% 2. FM 鉴频（使用更稳健的共轭相乘法）iq_data=iq_data-mean(iq_data);% 去除直流% 角度差分即瞬时频率fm_demod=angle(iq_data(2:end).*conj(iq_data(1:end-1)));baseband=fm_demod/(2*pi)*Fs;baseband=baseband/max(abs(baseband));% 归一化幅度，便于后续处理% 3. 设计滤波器lpFilt=designfilt('lowpassfir','FilterOrder',100,...'CutoffFrequency',15e3,'SampleRate',Fs);bpFilt_19k=designfilt('bandpassfir','FilterOrder',200,...'CutoffFrequency1',18.8e3,'CutoffFrequency2',19.2e3,...'SampleRate',Fs);bpFilt_38k=designfilt('bandpassfir','FilterOrder',200,...'CutoffFrequency1',37.8e3,'CutoffFrequency2',38.2e3,...'SampleRate',Fs);bpFilt_23_53=designfilt('bandpassfir','FilterOrder',200,...'CutoffFrequency1',23e3,'CutoffFrequency2',53e3,...'SampleRate',Fs);% 4. 提取各分量 (使用 filtfilt 进行零相位滤波，消除群延迟导致的相位失配)% 提取 L+RL_plus_R=filtfilt(lpFilt,baseband);% 提取 19kHz 导频pilot_19k=filtfilt(bpFilt_19k,baseband);% 恢复 38kHz 副载波（平方 + 带通滤波）subcarrier_38k=pilot_19k.*pilot_19k;subcarrier_38k=filtfilt(bpFilt_38k,subcarrier_38k)*2;% 乘2补偿幅度衰减% 提取并解调 L-R（23-53 kHz 带通 → 乘以副载波 → 低通）dsb_sc=filtfilt(bpFilt_23_53,baseband);L_minus_R=filtfilt(lpFilt,dsb_sc.*subcarrier_38k);% 5. 立体声矩阵解码L=(L_plus_R+L_minus_R)/2;R=(L_plus_R-L_minus_R)/2;% 6. 归一化右声道R=R/max(abs(R));% 7. 降采样到 48kHz 并播放（解决高采样率导致的声卡设备错误）Fs_play=48e3;R_play=resample(R,Fs_play,Fs);fprintf('准备播放右声道...\n');trysound(R_play,Fs_play);pause(length(R_play)/Fs_play+0.5);% 等待播放完毕catchMEfprintf('音频播放失败：%s\n',ME.message);end% 无论如何都保存一份 wav 方便试听audiowrite('R_audio.wav',R_play,Fs_play);fprintf('右声道音频已保存为 R_audio.wav (48kHz)\n');% 8. 自动解码摩斯电码（能量检测法）fprintf('\n--- 自动解码摩斯码 ---\n');% 发送端参数morse_symbol_len=floor(Fs*0.15);% 每个符号的采样点数（36000）% 对右声道信号进行 800 Hz 带通滤波，抑制带外噪声bp800=designfilt('bandpassiir','FilterOrder',4,...'HalfPowerFrequency1',750,'HalfPowerFrequency2',850,...'SampleRate',Fs);R_filt=filtfilt(bp800,R);% 零相位滤波% 计算短时能量（每符号长度一个窗口）num_symbols=12;% 已知摩斯码共 12 个符号energy=zeros(num_symbols,1);fork=1:num_symbols start_idx=(k-1)*morse_symbol_len+1;end_idx=min(k*morse_symbol_len,length(R_filt));segment=R_filt(start_idx:end_idx);energy(k)=mean(segment.^2);% 平均功率end% 动态阈值threshold=max(energy)*0.75;detected=double(energy>threshold);% 1: 有声, 0: 无声fprintf('检测到的摩斯码 ON/OFF 序列：');fprintf('%d ',detected);fprintf('\n');% 发送端原始模式（作为验证）expected_pattern=[101110111000];ifisequal(detected(:),expected_pattern(:))fprintf('摩斯码匹配！隐藏的 Flag 是：FLAG{FM_STEREO_R}\n');elsefprintf('模式不完全匹配，请人工收听 R_audio.wav 确认。\n');end% 绘制能量图用于观察figure;bar(energy);hold on;yline(threshold,'r--','Threshold');xlabel('Symbol index');ylabel('Energy');title('R 声道 800 Hz 能量分布');

七、 总结与预告

本题的真正难点不在于解调公式，而在于当实测频谱与理论描述严重不符时，能否系统性地诊断出差异来源。从采样率不足导致的混叠爆发，到三角噪声掩盖导频，再到高频时域波形的视觉错觉，每一个坑都可能让人误入歧途。正确的做法是：计算理论带宽验证采样率、用pwelch平滑频谱、用带通滤波做局部验证，让多个独立证据相互印证。
这套“频谱给方向、带通给证据”的方法论，在后续题目中会反复用到。下一题Q03 FM + RDS中，57 kHz 的 BPSK 副载波功率比 19 kHz 导频还低，单凭 FFT 几乎看不到。届时我们将引入锁相环（PLL / Costas Loop）——它不仅能从噪声中提取纯净相位，还能为 BPSK 的相干解调提供本地载波，是从“看谱”迈向“闭环跟踪”的关键一步。敬请期待。