其实在之前的在python下通过RTLSDR收听FM广播文章中, 已经找到了3种不同的FM解调方法, 只是着急得到结果, 所以没有详细研究, 现在已经搞定收音, 因此仔细研究一下解调方法的区别.

根据FM wikipedia的描述, 一般广播的调制范围为75kHz, 因此需要的采样率最低是150kHz, 那么我之前用的125kHz可以说是不恰当的, 由于是调频, 则波形应该是会受到削波.

其实对于已有IQ信号的情况而言, FM解调应当是比较容易的, 因为IQ信号天然的就具备的相位的信息, 把每个sample放在极坐标图上, 每个相邻sample之间角度差, 就是瞬时频率(关键词: instantaneous frequency), 把这个频率输出, 就得到了解调的信号. 那么这样就有了第一个方法 求反正切->求差->unwrap.

这里面有个问题, "求差"这个运算在连续时间的领域是微分, 关于在离散域做差分运算, 推荐参考Richard Lyons写的Novel Differentiator Handles Discrete Time-Domain Signals 或者买一本数字信号处理-(第三版) 莱昂斯 ISBN:9787121243677 看看第七章, 这本书第二版目前绝版, 而第三版是刚出的(2015年5月), 不知道翻译的怎么样.
"微分是很容易产生噪声", 看看模拟电路的微分器, 所以在离散运算中如果用最简单的做差(first-difference differentiator) 会带来两个问题:

  1. 数据不同步, 如果微分后的数据需要和原数据做计算, 则由于做差相当于一个(0.5,-0.5)两个tap的FIR滤波器, 那么它的group delay为0.5个sample(来自Lyons的书).
  2. 对高频信号进行了放大.

因此就有了书中介绍的一堆简单的FIR滤波器.

OK, 绕了一圈回到FM解调上, 可以看到上面这个方法需要加上一个unwrap, 这是由于角度不连续造成的, 在在python下实时显示rtlsdr波形与频谱中有截图. 那么如果换一种方法, 先求差, 然后再求角度, 这样是否就避免了这个问题呢? 在阅读librtlsdr中的rtl_fm程序中, 该程序就是采取的这个方法, 用当前sample乘后边sample的共轭, 然后再求atan2, 这样子每个角度都是从0度开始, 则不超过360度都不需要特殊处理.

还有一种方法是在Lyons的书中发现的, 位于13章的第22个, 在进行了一些导数的固定运算后得出了一个奇怪的式子.

FM解调公式

说实话, 我在这东西上面卡了好久, 推导过程是清晰明了的, 我在实际程序里也试用效果和上面两个差不多, 在网上翻了几天没啥发现, 我觉得自己研究一下, 于是我自造了几个数, 让这几个算法去解调.

    data = np.ones(16,np.complex)
    half= np.power(2,.5)/2
    data[0] = 1
    data[1] = half+half*1j
    data[2] = 1j
    data[3] = -half+half*1j
    data[4] = -1
    data[5] = -half-half*1j
    data[6] = -1j
    data[7] = half-half*1j
    data[8] = 1
    data[9] = half+half*1j
    data[10] = 1j
    data[11] = -half+half*1j
    data[12] = -3               #注意这里做的幅度不同的数
    data[13] = -half-half*1j
    data[14] = -5j              #同上
    data[15] = half-half*1j
    ...

为了省(tou)事(lan), 我造的是最简单的45度角, 绕着圆心转. 结果让我研究了两天.

#atan2->求差->unwrap    方案1
[ 0.          0.78539816  0.78539816  0.78539816  0.78539816  0.78539816
  0.78539816  0.78539816  0.78539816  0.78539816  0.78539816  0.78539816
  0.78539816  0.78539816  0.78539816  0.78539816]
#求角度差->atan2        方案2
[ 0.          0.78539816  0.78539816  0.78539816  0.78539816  0.78539816
  0.78539816  0.78539816  0.78539816  0.78539816  0.78539816  0.78539816
  0.78539816  0.78539816  0.78539816  0.78539816]
#专用公式+直接做差      方案3
[ 0.70710678  0.70710678  0.70710678  0.70710678  0.70710678  0.70710678
  0.70710678  0.70710678  0.70710678  0.70710678  0.70710678  0.23570226
  2.12132034  0.14142136  3.53553391]

可以看到前两个方案很完美的解出了所求的pi/4这个45度值(pi/4为0.78539816...). 专用公式为啥算成这个样子了? 而且还受幅度干扰严重, 前几个只所以能得到0.707这个貌似接近的值, 也是由 half= np.power(2,.5)/2 来提供的, 如果换成1, 就变成.

[ 0.          0.78539816  0.78539816  0.78539816  0.78539816  0.78539816
  0.78539816  0.78539816  0.78539816  0.78539816  0.78539816  0.78539816
  0.78539816  0.78539816  0.78539816  0.78539816]
[ 0.          0.78539816  0.78539816  0.78539816  0.78539816  0.78539816
  0.78539816  0.78539816  0.78539816  0.78539816  0.78539816  0.78539816
  0.78539816  0.78539816  0.78539816  0.78539816]
[ 0.5         1.          0.5         1.          0.5         1.          0.5
  1.          0.5         1.          0.5         0.33333333  1.5         0.2
  2.5       ]

完全不对了好吧, 于是我手算了一下, 这确实是算不出其他算法那样的结果, 一点儿可能性都没有, 不是1就是0, 怎么能算出来个pi/4 ? 于是我在困惑中怀疑上了那个微分运算, 也是在这样的情况下找到的那本书的第三版. 然而问题不是出在微分运算, 实验了多个并不完美的differentiator后, 依然得不到接近的结果. 于是只能google了.

Paul Solomon 在comp.dsp邮件组中多次关于关于他的FPGA FM收音机发问.

其他

说实话我看得不是很懂, 大概记得有人提到说是那个公式是一个小角度近似, 但我也没想明白为啥, 于是我决定试试. 首先需要做一个生成任意角度数据的.

    ticks = np.linspace(0,2*np.pi,10)   #角度由0至2*pi分成10等份.
    imag=np.cos(ticks)                  
    real=np.sin(ticks)*1j
    data = real+imag

首先试一下45度

#45度
ticks = np.linspace(0,np.pi*2,9)
[ 0.          0.78539816  0.78539816  0.78539816  0.78539816  0.78539816  0.78539816  0.78539816  0.78539816]
[ 0.          0.78539816  0.78539816  0.78539816  0.78539816  0.78539816  0.78539816  0.78539816  0.78539816]
[ 0.70710678  0.70710678  0.70710678  0.70710678  0.70710678  0.70710678  0.70710678  0.70710678]
#30度
ticks = np.linspace(0,np.pi*2,13)
[ 0.          0.52359878  0.52359878  0.52359878  0.52359878  0.52359878  0.52359878  0.52359878  0.52359878  0.52359878  0.52359878  0.52359878  0.52359878]
[ 0.          0.52359878  0.52359878  0.52359878  0.52359878  0.52359878  0.52359878  0.52359878  0.52359878  0.52359878  0.52359878  0.52359878  0.52359878]
[ 0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5]
#10度
ticks = np.linspace(0,np.pi*2,37)
[ 0.          0.17453293  0.17453293  0.17453293  0.17453293  0.17453293...  0.17453293]
[ 0.          0.17453293  0.17453293  0.17453293  0.17453293  0.17453293...  0.17453293]
[ 0.17364818  0.17364818  0.17364818  0.17364818  0.17364818  0.17364818...  0.17364818  0.17364818]
#90度
ticks = np.linspace(0,np.pi*2,5)
[ 0.          1.57079633  1.57079633  1.57079633  1.57079633]
[ 0.          1.57079633  1.57079633  1.57079633  1.57079633]
[ 1.  1.  1.  1.]

可以看到该公式随着角度的增加, 误差逐渐增大, 为啥就不清楚了. 也就是说, 该公式还是有局限性, 当采样率高时, 每个sample之间角度小, 就可以用, 如果角度大, 还是老老实实做atan吧~


Comments

comments powered by Disqus