利用USRP探索软件无线电(4)

2023-11-14 13:27| 来源: 网络整理| 查看: 265

21、引言

上回我们对IQ格式下的载波、AM调制和解调进行了展示。在此基础上，本篇继续研究另一种常用的模拟调制方式：FM。与AM恰好相反，FM调制方式只改变载波的相位，不改变载波的幅度。因此我们可以合理的推测，FM调制产生的复信号在复平面上应当只做切向运动，没有径向运动。

2、调频广播

调频广播是FM最开始最广泛的应用，也日常生活中常见的无线信号之一。除了日本之外，世界范围内的调频广播频段一般都在88~108MHz。调频广播有点类似移动通信的蜂窝结构，每个城市作为一个蜂窝小区。由于FM广播发射的功率是如此之强，导致经常能够收到相邻城市的信号。例如北京的FM广播频段看起来可能是这样的......

2.1 规范

大部分通信系统都是用规范(specification)来描述的，它描述了一个通信系统的业务需求、主要架构、各层的主要功能和接口，使用的关键技术和指标等。规范主要关注通信系统的架构和互操作的部分，一般不涉及具体的技术实现。它的作用类似建筑图纸，只描述一个通信系统应该是什么样子，至于如何达到，由施工单位，也就是通信设备厂商来完成。通过遵循预定义的、设计良好的接口，不同层次、厂商的通信设备无缝构成了一个完整的通信系统。

通信系统的规范一般由参与通信的各利益相关方共同制定。大部分的规范是公开的，但不是免费的，你需要花钱下载这个规范的pdf文档。自定义的通信规范一般是保密的，分析这种规范的行为被称为射频逆向工程(RF Reverse Engineering)。

我国现行的调频广播规范文件是《GB-T 4311-2000 米波调频广播技术规范》。

2.2 消息信号

FM广播虽然很常见，但是并不简单。为了传输立体声音频和其他信息，FM广播的消息信号使用了复用技术。下图显示了调频立体声广播的消息信号，图片来自参考[1]。

图1 立体声广播的消息信号

消息信号从左到右依次是：3-15kHz的“L+R”音频、19kHz的立体声导频，以38kHz为中心的DSB-SC“L-R”音频，以57kHz为中心的RDS信号，以67kHz或76kHz为中心的附加节目。立体声广播采取了一种“向下兼容”的方式，在保证了基础功能的前提下，为更多业务提供拓展支持。

3、分析FM信号 3.1 单音FM调制

本篇用到的原始信号在此下载。

上一篇提到了接收机频偏的影响，在接收18900kHz的信号时大约存在+10Hz的频偏，说明我的接收机晶振振荡偏慢，推算出偏差约为10/18900000=0.53ppm。我们将此值设置为GQRX的Freq.correction，这样GQRX会为某个具体的接收频率自动设置频偏。

与AM广播信号不同，FM是一种非线性调制。经过FM调制以后，已调信号中会产生消息信号中原本不存在的频率成分，这导致分析FM的调制规律比较复杂。大部分书中都用单音FM调制来分析FM调制的效果，例如参考[2]。

上一篇我们利用AM广播的静音间隙得到了载波的波形。对于FM立体声广播，消息信号中始终存在19kHz的单音，因此在静音时我们会得到单音(19kHz)FM调制的波形。有一些电台每周二下午进行设备检修，此时会发射静音信号。

利用下面的代码画出单音FM调制的波形：

%read samples fid=fopen('gqrx_20211109_085354_104400065_250000_fc.raw'); x = fread(fid,'float'); x = x(1:2:end) + 1i*x(2:2:end); lower=10000; upper=10250; rp=real(x(lower:upper)); ip=imag(x(lower:upper)); sa=(1:upper-lower+1); % plot3(rp,ip,sa,'Marker','.'); lim=0.005; axis([-lim,lim,-lim,lim,1,length(sa)]); ax = gca; ax.XAxisLocation = 'origin'; ax.YAxisLocation = 'origin'; xlabel('real'); ylabel('imaginary'); zlabel('samples'); grid on; view(3); waveform = animatedline('Color','b','Marker','.'); motion = animatedline('Color','r','Marker','.'); for k = 1:length(sa) %clearpoints(motion); addpoints(motion,rp(k),ip(k),1); addpoints(waveform,rp(k),ip(k),sa(k)); drawnow end view(0,90);

从下图可以看到用正弦消息信号对复载波进行FM调制的结果。我们选取了250个样本点，其持续时间为250/250000=1ms。对于19kHz的正弦信号，1毫秒内会经历19个完整周期。当正弦信号的幅度为0时，对应已调信号拐弯的地方。当正弦信号到达波峰或者波谷，已调信号旋转运动到周期的中间，此时角速度的绝对值最大。通过计算相邻样本点的相位差，再乘采样率250000，就能得到该时间内FM调制的频偏。

3.2 FM解调

FM解调即求解已调信号的瞬时频率的过程。直观的想法是先用arctan(imag/real)求出每个样本点的相位，然后对相邻样本点的相位做差分，即可得到频率。但是用atan或atan2求反正切会遇到相位展开的问题。所以在实践中都采用另外一种算法，即先求 $tan(\theta _{n}-\theta_{n-1})$ ，再求反正切。前者可以用当前样本点和上一个样本点的共轭相乘得到，即：

$[cos(\theta_{n})+jsin(\theta_{n})]*[cos(\theta_{n-1})-jsin(\theta_{n-1})] =$ $cos(\theta_{n}-\theta_{n-1})+jsin(\theta_{n}-\theta_{n-1})$

这里还有不用求反正切的FM解调算法。

106.1MHz中国之声的FM解调代码如下：

fid=fopen('gqrx_20211109_090124_106100063_250000_fc.raw'); x = fread(fid,'float'); x = x(1:2:end) + 1i*x(2:2:end); fs=250000; %工具箱FM解调 fd=75; DEMOD = comm.FMDemodulator('SampleRate',fs,'FrequencyDeviation',fd); z0 = step(DEMOD,x); %手工解调 delay=[x(2:length(x));0]; c=delay.*conj(x); z1=angle(c); %选取一个片段画频谱 n=8192; origin=1000; window=hanning(n); s0=z1(origin:origin+n-1); s0=s0.*window; y=fft(s0); y0=fftshift(y); f0 = (-n/2:n/2-1)*(fs/n); power=abs(y0).^2/n; logpower=10*(log10(power)); plot(f0,logpower); xlim([-fs/2,fs/2]); %写入文件 audiowrite('fm.wav',z1,fs,'BitsPerSample',32); %工具箱画频谱 SA = dsp.SpectrumAnalyzer('SampleRate',fs,'ShowLegend',true); step(SA,[z0 z1]);

两种方法画出的功率谱图：

用Adobe Audition打开生成的fm.wav文件，按alt+z查看频谱：

用500kHz带宽查看106.1MHz的频谱如下，可以说这个广播把能用的技术都用上了：立体声、RDS、SCA，以及CDR(Chinese Digital Radio，HD Radio中国版)。

3.3 立体声解调

完成第一层的FM解调之后，继续进行立体声解调，大体思路是分别滤波，得到L+R、19kHz导频，L-R，利用导频进行DSB-SC解调恢复L-R，在通过加减运算得到L和R。

【本文地址】

公司简介

联系我们