信息处理 - 使用相位导数的 IQ 信号 FM 判别 - 吾爱随笔录

我正在尝试解调以 2.4Msps 从 RTL-SDR 记录为 IQ 样本的 FM 信号，其中 RTL-SDR 以 FM 电台的频率为中心。

所以输入信号是 2.4MSps 的 8 位 IQ 基带。

数据可在此处获得。

输入信号的 FFT 看起来很好，与我在 SDR-Sharp 中看到的相匹配，所以我认为数据很好。

我的理解是，我应该能够通过以下方式区分 Octave 中的 FM：

消除直流偏移
将其低通滤波至 240KHz 并抽取
计算相位分量的离散导数

FFT 在低通滤波和抽取之后看起来不错。
我的假设是应用FIR 低通滤波器不应该扭曲相位（不确定我是否正确）

然后，我通过将相位信号与 [1,-1] 的离散导数的脉冲响应进行卷积来计算信号相位的离散导数。

我期望衍生信号的频谱将包含 0-15KHz 的复合 FM 信号，由单声道 L+R 分量、19KHz 导频音等占据。
但是我在相位导数的 FFT 中得到白噪声（见下文）。

我的八度代码如下

% 2.4MSps IQ baseband signal of an FM station from RTL-SDR recorded in SDR#
  samples = audioread("SDRSharp_20191229_133205Z_102102000Hz_IQ.wav", "native");
  csamples = arrayfun(@complex,samples(:,1),samples(:,2));

% 1) - remove DC offset
  csamples_ac = csamples - mean(csamples);   

% 2) - low pass filter and decimate to get complex 240KSps (240KHz Widband FM)
  csamples_filtered = decimate(csamples_ac, 10, "fir");

% 3) - descrete derivative of phase
  filtered_derivative = conv(angle(csamples_filtered), [1,-1]);

  plot(abs(fftshift(fft(filtered_derivative, 32768))));
  axis("tight");

编辑：感谢 hotpaw2 和 Dan Boschen，我现在知道我的错误是什么，并且有了一个有效的 Octave 脚本。主要的错误是没有展开相位，同样当通过卷积与 [1,-1] 进行区分时，需要丢弃结果中的第一个和最后一个样本，因为卷积的可能伪影会在序列的末端给出巨大的值，第三个错误是在使用 Octave 的音频播放器功能之前没有转换为 int16 并且没有去除声音样本中的 DC。

如果有人需要这样的脚本，下面是一个工作版本。
（请注意，基于相位的鉴别会受到低通滤波 1/f 滚降效应的影响，有关更多详细信息，请参阅此答案）

pkg load signal;

#reading SDR Sharp recording (2.4Msps,RAW,8bit,baseband)
samples = audioread("SDRSharp_20191229_183956Z_99500000Hz_IQ.wav", "native");
#one channel real, other channel - imaginary components
unshifted_signal = double(samples(:,1)) .+ double(samples(:,2))*j;

#removing DC, this is important, 
unshifted_signal = unshifted_signal - mean(unshifted_signal);

#a 1550 Hz frequency shift to move the real carrier center frequency as close to 0 as possible (de-trending the phase)
#the frequency shift is not absolutely necessary, small diviations don't affect the demodulation
#this would be station specific, selected to minimized the slope on plot(unwrp_phase)
freq_shift = exp(-j*1550*2*pi*[1:length(unshifted_signal)]/2400000)';
#freq_shift = 1;
signal = unshifted_signal.*freq_shift;

#low-pass filter and decimate down to 240KHz
rcv_240 = decimate(signal, 10);
#calculating the phase, and unwrapping it (unwrapping is important)
unwrp_phase = unwrap(angle(rcv_240));

#doing the actual demodulation (discrimination)
#not including the ends because of possible convolution artifacts
phase_drv = conv(unwrp_phase, [1,-1])(2:end-1);

#extract L+R submodulated at 0-15KHz by low-pass filtering and decimation
monoLplusR_unnormalized = decimate(phase_drv, 8);

#normalize the range for 16bits for audioplayer, it is important to have no DC here (removed earlier)
monoLplusR = int16(monoLplusR_unnormalized/max(abs(monoLplusR_unnormalized))*32766);

#play at 30Ksps (0-15KHz)
player = audioplayer(monoLplusR, 30000, 16);
play(player);