说明
学习数字信号处理算法时整理的学习笔记。同系列文章目录可见 《DSP 学习之路》目录,代码已上传到 Github - ModulationAndDemodulation。本篇介绍 FM 调频信号的调制与解调,内附全套 MATLAB 代码。
文章目录
- 说明
- 1. FM 调制算法
- 2. FM 解调算法
- 3. FM 仿真(MATLAB Communications Toolbox)
- 参考资料
- 附录代码
- 附.1 文件 mod_fm.m
- 附.2 文件 main_modFM_example1.m
- 附.3 文件 main_modFM_example2.m
- 附.4 文件 lpf_filter.m
- 附.5 文件 demod_fm_method1.m
- 附.6 文件 main_demodFM_example1.m
- 附.7 文件 demod_fm_method2.m
- 附.8 文件 main_demodFM_example2.m
- 附.9 文件 demod_fm_method3.m
- 附.10 文件 main_demodFM_example3.m
- 附.11 文件 demod_fm_method4.m
- 附.12 文件 main_demodFM_example4.m
- 附.13 文件 main_CommFM_example.m
1. FM 调制算法
1.1 FM 信号描述
用调制信号去控制载波的瞬时频率,使其按照调制信号的规律变化,当调制信号是模拟信号时,这个过程就被称为调频(FM)。FM 信号的时域表达式为:
s
F
M
(
t
)
=
A
c
o
s
[
ω
c
t
+
K
f
∫
0
t
m
(
τ
)
d
τ
]
(1)
s_{FM}(t)=Acos{\left[\omega_ct+{K_f}\int_0^{t}{m(\tau)}d{\tau}\right]} \tag{1}
sFM(t)=Acos[ωct+Kf∫0tm(τ)dτ](1)
式中:
A
A
A 为载波恒定振幅,
K
f
K_f
Kf 为调频灵敏度(单位
r
a
d
/
(
s
⋅
V
)
rad/(s{\cdot}V)
rad/(s⋅V)),
m
(
t
)
m(t)
m(t) 是调制信号(携带要发出去的信息),
c
o
s
ω
c
t
cos{\omega_ct}
cosωct 是载波,
ω
c
\omega_c
ωc 是载波角频率,与载波频率
f
c
f_c
fc 之间的关系为
ω
c
=
2
π
f
c
\omega_c=2{\pi}f_c
ωc=2πfc。由式
(
1
)
(1)
(1) 可得 FM 信号相对于载频
ω
c
{\omega_c}
ωc 的瞬时频偏为:
d
[
K
f
∫
0
t
m
(
τ
)
d
τ
]
d
t
=
K
f
m
(
t
)
(2)
\frac{d{\left[{K_f}\int_0^{t}{m(\tau)}d{\tau}\right]}}{dt}={K_f}m(t) \tag{2}
dtd[Kf∫0tm(τ)dτ]=Kfm(t)(2)
由
(
1
)
(1)
(1) 式可知 FM 信号相对于载波相位
ω
c
t
{\omega}_ct
ωct 的瞬时相位偏移随
m
(
t
)
m(t)
m(t) 的积分呈线性变化,由
(
2
)
(2)
(2) 式可知 FM 信号相对于载频
ω
c
{\omega_c}
ωc 的瞬时频率偏移随
m
(
t
)
m(t)
m(t) 成线性变化,比例系数都为
K
f
K_f
Kf。有时候会用
k
f
k_f
kf 表示调频灵敏度(单位
H
z
/
V
Hz/V
Hz/V),它与
K
f
K_f
Kf 之间的关系为
K
f
=
2
π
k
f
{K_f}=2{\pi}{k_f}
Kf=2πkf,需注意这个转换关系。FM 的调频指数(调制指数)
β
{\beta}
β 定义如下,其中
W
W
W 是基带信号(调制信号)
m
(
t
)
m(t)
m(t) 的带宽(或最高频率):
β
=
k
f
∣
m
(
t
)
∣
m
a
x
W
(3)
{\beta}=\frac{{k_f}{\lvert}{m(t)}{\rvert}_{max}}{W} \tag{3}
β=Wkf∣m(t)∣max(3)
若
m
(
t
)
m(t)
m(t) 为单一频率的正弦波(即
m
(
t
)
=
A
m
c
o
s
(
2
π
f
m
t
)
m(t)={A_m}cos({2{\pi}{f_m}t})
m(t)=Amcos(2πfmt) 时,此时
W
=
f
m
W={f_m}
W=fm),则调制指数的表达式如下,调制指数的值正好与最大相位偏移相同,其中
Δ
f
=
β
f
m
{\Delta}f={\beta}{f_m}
Δf=βfm 表示最大频偏。
β
=
K
f
A
m
ω
m
=
k
f
A
m
f
m
=
Δ
f
f
m
=
[
K
f
∫
0
t
m
(
τ
)
d
τ
]
m
a
x
(4)
{\beta}=\frac{{K_f}{A_m}}{\omega_m}=\frac{{k_f}{A_m}}{f_m}=\frac{{\Delta}f}{f_m}=\left[{K_f}\int_0^{t}{m(\tau)}d{\tau}\right]_{max} \tag{4}
β=ωmKfAm=fmkfAm=fmΔf=[Kf∫0tm(τ)dτ]max(4)
FM 信号的频域表达式比较复杂,下面分成窄带调频与宽带调频两种情况来简化讨论。
(1)窄带调频
一般将由
m
(
t
)
m(t)
m(t) 引起的最大瞬时相位偏移远小于
3
0
∘
30^{\circ}
30∘ 的情况称为窄带 FM,即满足以下条件时,FM 信号的频谱宽度比较窄,称为窄带调频(NBFM)。窄带调频占据的带宽较窄,传输数据量有限,主要应用于无线语音的传输(比如无线对讲机)。
∣
K
f
∫
0
t
m
(
τ
)
d
τ
∣
≪
π
6
(5)
{\lvert}{K_f}\int_0^{t}{m(\tau)}d{\tau}{\rvert}\ll\frac{\pi}{6} \tag{5}
∣Kf∫0tm(τ)dτ∣≪6π(5)
此时,式
(
1
)
(1)
(1) 可以近似为:
s
N
B
F
M
(
t
)
=
A
c
o
s
[
ω
c
t
+
K
f
∫
0
t
m
(
τ
)
d
τ
]
=
A
c
o
s
(
ω
c
t
)
c
o
s
[
K
f
∫
0
t
m
(
τ
)
d
τ
]
−
A
s
i
n
(
ω
c
t
)
s
i
n
[
K
f
∫
0
t
m
(
τ
)
d
τ
]
≈
A
c
o
s
(
ω
c
t
)
⋅
1
−
A
s
i
n
(
ω
c
t
)
⋅
K
f
∫
0
t
m
(
τ
)
d
τ
=
A
c
o
s
(
ω
c
t
)
−
[
A
K
f
∫
0
t
m
(
τ
)
d
τ
]
s
i
n
(
ω
c
t
)
(6)
\begin{aligned} s_{NBFM}(t)&=Acos{\left[\omega_ct+{K_f}\int_0^{t}{m(\tau)}d{\tau}\right]}\\[1em] &=Acos{(\omega_ct)}cos{\left[{K_f}\int_0^{t}{m(\tau)}d{\tau}\right]}-Asin{(\omega_ct)}sin{\left[{K_f}\int_0^{t}{m(\tau)}d{\tau}\right]}\\[1em] &\approx{Acos}{(\omega_ct)}\cdot{1}-Asin{(\omega_ct)}\cdot{{K_f}\int_0^{t}{m(\tau)}d{\tau}}\\[1em] &={Acos}{(\omega_ct)}-\left[A{K_f}\int_0^{t}{m(\tau)}d{\tau}\right]sin{(\omega_ct)} \end{aligned} \tag{6}
sNBFM(t)=Acos[ωct+Kf∫0tm(τ)dτ]=Acos(ωct)cos[Kf∫0tm(τ)dτ]−Asin(ωct)sin[Kf∫0tm(τ)dτ]≈Acos(ωct)⋅1−Asin(ωct)⋅Kf∫0tm(τ)dτ=Acos(ωct)−[AKf∫0tm(τ)dτ]sin(ωct)(6)
对其做傅里叶变换,得到窄带调频(NBFM)信号的频谱(幅度谱)表达式:
S
N
B
F
M
(
ω
)
=
π
A
[
δ
(
ω
+
ω
c
)
+
δ
(
ω
−
ω
c
)
]
+
A
K
f
2
[
M
(
ω
−
ω
c
)
ω
−
ω
c
−
M
(
ω
+
ω
c
)
ω
+
ω
c
]
(7)
S_{NBFM}(\omega)={\pi}A\left[\delta(\omega+\omega_c)+\delta(\omega-\omega_c)\right]+\frac{AK_f}{2}\left[\frac{M(\omega-\omega_c)}{\omega-\omega_c}-\frac{M(\omega+\omega_c)}{\omega+\omega_c}\right] \tag{7}
SNBFM(ω)=πA[δ(ω+ωc)+δ(ω−ωc)]+2AKf[ω−ωcM(ω−ωc)−ω+ωcM(ω+ωc)](7)
式中,
M
(
ω
)
M(\omega)
M(ω) 是调制信号
m
(
t
)
m(t)
m(t) 的频谱。与 AM 信号不同的是,NBFM 信号的两个边频分别乘了因式
1
/
(
ω
−
ω
c
)
1/(\omega-\omega_c)
1/(ω−ωc) 和
1
/
(
ω
+
ω
c
)
1/(\omega+\omega_c)
1/(ω+ωc),由于因式是频率的函数,所以这种加权是频率加权,加权的结果引起调制信号频谱的失真,此外, NBFM 的一个边带和 AM 反相。
(2)宽带调频
当不满足
(
5
)
(5)
(5) 式所表示的条件时,FM 信号被称为宽带调频(WBFM)。宽带调频所占用的频带宽度比较宽,传输数据量大,主要应用于调频立体声广播。WBFM 的时域表达式无法进一步简化,首先考虑单音调制的情况,然后把分析的结论推广到多音调制,当
m
(
t
)
=
A
m
c
o
s
(
ω
m
t
)
m(t)={A_m}cos({{\omega_m}t})
m(t)=Amcos(ωmt) 时,带入式
(
1
)
(1)
(1) 展开可得:
s
W
B
F
M
(
t
)
=
A
c
o
s
[
ω
c
t
+
K
f
∫
0
t
A
m
c
o
s
(
ω
m
τ
)
d
τ
]
=
A
c
o
s
[
ω
c
t
+
K
f
A
m
ω
m
s
i
n
(
ω
m
t
)
]
=
A
c
o
s
[
ω
c
t
+
β
s
i
n
(
ω
m
t
)
]
=
A
∑
n
=
−
∞
∞
J
n
(
β
)
c
o
s
[
(
ω
c
+
n
ω
m
)
t
]
(8)
\begin{aligned} s_{WBFM}(t)&=Acos{\left[\omega_ct+{K_f}\int_0^{t}{{A_m}cos({{\omega_m}{\tau}})}d{\tau}\right]}\\[1em] &=Acos{\left[\omega_ct+\frac{{K_f}{A_m}}{\omega_m}sin({{\omega_m}t})\right]}\\[1em] &=Acos{\left[\omega_ct+{\beta}sin({{\omega_m}t})\right]}\\[1em] &=A\sum\limits_{n=-\infty}^{\infty}{{J_n}(\beta)cos{\left[(\omega_c+n\omega_m)t\right]}} \end{aligned} \tag{8}
sWBFM(t)=Acos[ωct+Kf∫0tAmcos(ωmτ)dτ]=Acos[ωct+ωmKfAmsin(ωmt)]=Acos[ωct+βsin(ωmt)]=An=−∞∑∞Jn(β)cos[(ωc+nωm)t](8)
式中
J
n
(
β
)
J_n(\beta)
Jn(β) 为第一类
n
n
n 阶贝塞尔函数,它是调频指数
β
\beta
β 的函数。对其做傅里叶变换,得到单音调制时宽带调频(WBFM)信号的频谱(幅度谱)表达式:
S
W
B
F
M
(
ω
)
=
π
A
∑
n
=
−
∞
∞
J
n
(
β
)
[
δ
(
ω
−
ω
c
−
n
ω
m
)
+
δ
(
ω
+
ω
c
+
n
ω
m
)
]
(9)
S_{WBFM}(\omega)={\pi}A\sum\limits_{n=-\infty}^{\infty}{{J_n}(\beta)\left[\delta(\omega-\omega_c-n\omega_m)+\delta(\omega+\omega_c+n\omega_m)\right]} \tag{9}
SWBFM(ω)=πAn=−∞∑∞Jn(β)[δ(ω−ωc−nωm)+δ(ω+ωc+nωm)](9)
由式
(
8
)
(8)
(8) 和式
(
9
)
(9)
(9) 可见,WBFM 调频信号的频谱由载波分量
ω
c
\omega_c
ωc 和无数边频
ω
c
±
n
ω
m
{\omega_c}\pm{n\omega_m}
ωc±nωm 组成。当
n
=
0
n=0
n=0 时是载波分量
ω
c
\omega_c
ωc,其幅度为
A
J
0
(
β
)
AJ_0(\beta)
AJ0(β),当
n
≠
0
n\not=0
n=0 时就是对称分布在载频两侧的边频分量
ω
c
±
n
ω
m
{\omega_c}\pm{n\omega_m}
ωc±nωm,其幅度为
A
J
n
(
β
)
AJ_n(\beta)
AJn(β),相邻边频之间的间隔为
ω
m
\omega_m
ωm。且当
n
n
n 为奇数时,上下边频极性相反,当
n
n
n 为偶数时极性相同。由此可见,FM 信号的频谱不再是调制信号频谱的线性搬移,而是一种非线性过程。
1.2 FM 信号的带宽与功率分配
宽带调频信号的频谱包含无穷多个频率分量,因此理论上调频信号的频带宽度为无限宽。但是,实际上边频幅度
J
n
(
β
)
J_n(\beta)
Jn(β) 随着
n
n
n 的增大而逐渐减小,因此只要取适当的
n
n
n 值使边频分量小到可以忽略的程度,调频信号可以近似认为具有有限频谱。通常采用的原则是:信号的频带宽度应包括幅度大于未调载波的
10
%
10\%
10% 以上的边频分量,即
∣
J
n
(
β
)
∣
≥
0.1
\lvert{J_n(\beta)}\rvert\geq0.1
∣Jn(β)∣≥0.1。根据经验,当
β
≥
1
\beta\geq1
β≥1 时,取边频数
n
=
β
+
1
n=\beta+1
n=β+1 即可,因为
n
>
β
+
1
n>\beta+1
n>β+1 以上的边频幅度
J
n
(
β
)
J_n(\beta)
Jn(β) 均小于
0.1
0.1
0.1,相应产生的功率均在总功率的
2
%
2\%
2% 以下,可以忽略不计,根据这条经验规则,调频波的有效带宽为:
B
F
M
=
2
(
β
+
1
)
W
(10)
B_{FM}=2(\beta+1)W \tag{10}
BFM=2(β+1)W(10)
这就是广泛用于计算调频信号带宽的卡森(Carson)公式。式中
β
\beta
β 为调频指数,
W
W
W 为基带信号(调制信号)
m
(
t
)
m(t)
m(t) 的带宽(或最高频率)。对于单音调制信号(即
m
(
t
)
=
A
m
c
o
s
(
2
π
f
m
t
)
m(t)={A_m}cos({2{\pi}{f_m}t})
m(t)=Amcos(2πfmt) 时,此时
W
=
f
m
W={f_m}
W=fm),将
(
4
)
(4)
(4) 式带入
(
10
)
(10)
(10) 式,可得:
B
F
M
=
2
(
β
+
1
)
f
m
=
2
(
Δ
f
+
f
m
)
(11)
B_{FM}=2(\beta+1)f_m=2({\Delta}f+f_m) \tag{11}
BFM=2(β+1)fm=2(Δf+fm)(11)
- 当 β ≪ 1 \beta \ll 1 β≪1 时(对应窄带调频),带宽可近似估计为 B N B F M ≈ 2 f m B_{NBFM}\approx2f_m BNBFM≈2fm,这时,带宽由第一对边频分量决定,带宽只随调制频率 f m f_m fm 变化,而与最大频偏 Δ f {\Delta}f Δf 无关。
- 当 β ≫ 1 \beta \gg 1 β≫1 时(对应宽带调频),带宽可近似估计为 B W B F M ≈ 2 Δ f B_{WBFM}\approx2{\Delta}f BWBFM≈2Δf,这时,带宽由最大频偏 Δ f {\Delta}f Δf 决定,而与调制频率 f m f_m fm 无关。
利用帕塞瓦尔定理以及贝塞尔函数性质,不难求得 FM 信号功率
P
F
M
P_{FM}
PFM 与未调载波功率
P
c
P_c
Pc 的关系如下:
P
F
M
=
s
W
B
F
M
2
(
t
)
‾
=
A
2
2
∑
n
=
−
∞
∞
J
n
2
(
β
)
=
A
2
2
=
P
c
(12)
P_{FM}=\overline{s_{WBFM}^2(t)}=\frac{A^2}{2}\sum\limits_{n=-\infty}^{\infty}{{J_n^2}(\beta)}=\frac{A^2}{2}=P_c \tag{12}
PFM=sWBFM2(t)=2A2n=−∞∑∞Jn2(β)=2A2=Pc(12)
上式说明,调频信号的平均功率
P
F
M
P_{FM}
PFM 等于未调载波的平均功率
P
c
P_c
Pc,即调制后总的功率不变,只是将原来载波功率中的一部分分配给每个边频分量,所以,调制过程只是进行功率的重新分配,而分配的原则与调频指数
β
\beta
β 有关。
1.3 FM 信号的调制方法
FM 信号的调制方法有 3 种,一种是直接调频法,一种是间接调频法,第三种是正交调制法。
(1)直接调频法
直接调频法是用调制信号 m ( t ) m(t) m(t) 直接控制高频振荡器,让回路元件的参数发生改变,使其输出频率按调制信号的规律线性地变化,常用的元件是变容二极管。直接调频法的主要优点是在实现线性调频的要求下,可以获得较大的频偏,且实现电路简单;主要缺点是频率稳定度不高,往往需要采用自动频率控制系统来稳定中心频率。
(2)间接调频法
间接调频法也被称为倍频法、阿姆斯特朗法。该方法先将调制信号积分,然后对载波进行调相,即可产生一个 NBFM 信号,再经 n n n 次倍频器得到 WBFM 信号,间接调频法的优点是频率稳定性好,缺点是需要多次倍频和混频,电路较复杂。如下图所示,根据式 ( 6 ) (6) (6) 窄带调频的时域表达式可知,虚线框内所示部分可以用来获得 NBFM 信号。
上图中的倍频器是一个非线性器件,以理想平方律器件为例,其输出
x
o
(
t
)
x_o(t)
xo(t) 与输入
x
i
(
t
)
x_i(t)
xi(t) 的关系为
x
o
(
t
)
=
a
x
i
2
(
t
)
{x_o}(t)=a{x_i^2}(t)
xo(t)=axi2(t),其中
a
a
a 为常数,将式
(
6
)
(6)
(6) 带入,得到:
x
o
(
t
)
=
a
x
2
(
t
)
=
1
2
a
A
2
{
1
+
c
o
s
[
2
ω
c
t
+
2
K
f
∫
0
t
m
(
τ
)
d
τ
]
}
(13)
{x_o}(t)=ax^2(t)=\frac{1}{2}aA^2\left\{ 1+cos{\left[2\omega_ct+2{K_f}\int_0^{t}{m(\tau)}d{\tau}\right]}\right\} \tag{13}
xo(t)=ax2(t)=21aA2{1+cos[2ωct+2Kf∫0tm(τ)dτ]}(13)
接着将
x
o
(
t
)
x_o(t)
xo(t) 经过一个带通滤波器,就可以将其中的直流分量滤除,从而得到一个新的 FM 信号
x
o
′
(
t
)
x'_o(t)
xo′(t),如下:
x
o
′
(
t
)
=
1
2
a
A
2
c
o
s
[
2
ω
c
t
+
2
K
f
∫
0
t
m
(
τ
)
d
τ
]
(14)
x'_o(t)=\frac{1}{2}aA^2cos{\left[2\omega_ct+2{K_f}\int_0^{t}{m(\tau)}d{\tau}\right]} \tag{14}
xo′(t)=21aA2cos[2ωct+2Kf∫0tm(τ)dτ](14)
这个信号的载频和相位偏移均增为原来的
2
2
2 倍,由于相位偏移增为
2
2
2 倍,因而调频指数也必然增为
2
2
2 倍,同理,经
n
n
n 次倍频后可以使调频信号的载频和调频指数增为
n
n
n 倍。增大调制指数时,一般需要保持载频不变或者不增大太多,这个时候可以接着使用一个混频器配一个带通滤波器来进行下变频,混频器只改变载频而不影响频偏,具体实现方法可参考《通信原理》相关资料。
(3)正交调制法
将
(
1
)
(1)
(1) 式进行三角展开,可以得到:
s
F
M
(
t
)
=
A
c
o
s
[
ω
c
t
+
K
f
∫
0
t
m
(
τ
)
d
τ
]
=
A
c
o
s
(
ω
c
t
)
c
o
s
[
K
f
∫
0
t
m
(
τ
)
d
τ
]
−
A
s
i
n
(
ω
c
t
)
s
i
n
[
K
f
∫
0
t
m
(
τ
)
d
τ
]
(15)
\begin{aligned} s_{FM}(t)&=Acos{\left[\omega_ct+{K_f}\int_0^{t}{m(\tau)}d{\tau}\right]}\\[1em] &=Acos{(\omega_ct)}cos{\left[{K_f}\int_0^{t}{m(\tau)}d{\tau}\right]}-Asin{(\omega_ct)}sin{\left[{K_f}\int_0^{t}{m(\tau)}d{\tau}\right]}\\[1em] \end{aligned} \tag{15}
sFM(t)=Acos[ωct+Kf∫0tm(τ)dτ]=Acos(ωct)cos[Kf∫0tm(τ)dτ]−Asin(ωct)sin[Kf∫0tm(τ)dτ](15)
正交调制流程如下:
- 对调制信号 m ( t ) m(t) m(t) 进行积分,得到 Φ = K f ∫ 0 t m ( τ ) d τ \Phi=K_f\int_0^{t}{m(\tau)}d{\tau} Φ=Kf∫0tm(τ)dτ。
- 对积分后的信号分别取余弦和正弦,得到 I I I 路数据 I ( t ) = c o s ( Φ ) I(t)=cos(\Phi) I(t)=cos(Φ) 与 Q Q Q 路数据 Q ( t ) = s i n ( Φ ) Q(t)=sin(\Phi) Q(t)=sin(Φ)。
- 分别乘以载波 A c o s ( ω c t ) Acos(\omega_ct) Acos(ωct) 与 − A s i n ( ω c t ) -Asin(\omega_ct) −Asin(ωct) 后相加,得到 FM 信号 s F M ( t ) = I ( t ) A c o s ( ω c t ) − Q ( t ) A s i n ( ω c t ) s_{FM}(t)=I(t)Acos(\omega_ct)-Q(t)Asin(\omega_ct) sFM(t)=I(t)Acos(ωct)−Q(t)Asin(ωct)。
其中第三步也可以先将 I ( t ) I(t) I(t) 和 Q ( t ) Q(t) Q(t) 组成一个复信号 Z ( t ) = I ( t ) + i Q ( t ) Z(t)=I(t)+iQ(t) Z(t)=I(t)+iQ(t),然后乘以复载波 e x p ( i ω c t ) exp(i\omega_ct) exp(iωct),最后取实部,即 s F M ( t ) = R e a l [ Z ( t ) e x p ( i ω c t ) ] s_{FM}(t)=Real\left[Z(t)exp(i\omega_ct)\right] sFM(t)=Real[Z(t)exp(iωct)],两种方法是等价的。
1.4 窄带 FM 信号示例
上一小节中介绍的调制方法都与硬件实现相关,接下来使用 MATLAB 软件来直接生成 NBFM 信号。调制信号
m
(
t
)
m(t)
m(t) 可以是确知信号,也可以是随机信号。当
m
(
t
)
m(t)
m(t) 是确知信号时,不妨假设
m
(
t
)
m(t)
m(t) 的时域表达式如下,对应的基带带宽
W
=
f
m
W=f_m
W=fm:
m
(
t
)
=
s
i
n
(
2
π
f
m
t
)
+
c
o
s
(
π
f
m
t
)
(16)
m(t) = sin(2{\pi}{f_m}t)+cos({\pi}{f_m}t) \tag{16}
m(t)=sin(2πfmt)+cos(πfmt)(16)
各调制参数取值:
A
=
1
A=1
A=1,
f
m
=
2500
H
z
f_m=2500Hz
fm=2500Hz,
β
=
1
0
−
6
{\beta}=10^{-6}
β=10−6,
f
c
=
20000
H
z
f_c=20000Hz
fc=20000Hz。信号采样率
f
s
=
8
f
c
f_s=8{f_c}
fs=8fc,仿真总时长为
2
s
2s
2s。FM 调制效果如下图所示(为了美观,时域只显示前 500 个点),调制信号
m
(
t
)
m(t)
m(t) 双边幅度谱有四根离散谱线(
±
2500
H
z
{\pm}2500Hz
±2500Hz、
±
1250
H
z
{\pm}1250Hz
±1250Hz),调制信号
m
(
t
)
m(t)
m(t) 积分后的双边幅度谱有五根离散谱线(
0
H
z
0Hz
0Hz、
±
2500
H
z
{\pm}2500Hz
±2500Hz、
±
1250
H
z
{\pm}1250Hz
±1250Hz),NBFM 调频信号
s
(
t
)
s(t)
s(t) 双边幅度谱有十根离散谱线(
±
22500
H
z
{\pm}22500Hz
±22500Hz、
±
21250
H
z
{\pm}21250Hz
±21250Hz、
±
20000
H
z
{\pm}20000Hz
±20000Hz、
±
18750
H
z
{\pm}18750Hz
±18750Hz、
±
17500
H
z
{\pm}17500Hz
±17500Hz)。NBFM 调频信号
s
(
t
)
s(t)
s(t) 的带宽满足公式
B
N
B
F
M
≈
2
f
m
=
5000
H
z
B_{NBFM}\approx2f_m=5000Hz
BNBFM≈2fm=5000Hz,代码详见附录 main_modFM_example1.m
与 mod_fm.m
,设置 beta = 1e-6
即可。
当
m
(
t
)
m(t)
m(t) 是随机信号时,不妨假设基带信号带宽为
W
=
f
H
=
3000
H
z
W={f_H}=3000Hz
W=fH=3000Hz,各调制参数取值:
A
=
1
A=1
A=1,
β
=
1
0
−
6
{\beta}=10^{-6}
β=10−6,
f
c
=
20000
H
z
f_c=20000Hz
fc=20000Hz。信号采样率
f
s
=
8
f
c
f_s=8{f_c}
fs=8fc,仿真总时长为
2
s
2s
2s。FM 调制效果如下图所示(为了美观,时域只显示前 500 个点),调制信号
m
(
t
)
m(t)
m(t) 双边幅度谱中间谱峰的范围约为
−
3000
H
z
∼
3000
H
z
-3000Hz{\sim}3000Hz
−3000Hz∼3000Hz,调制信号
m
(
t
)
m(t)
m(t) 积分后的双边幅度谱中间谱峰的范围依然约为
−
3000
H
z
∼
3000
H
z
-3000Hz{\sim}3000Hz
−3000Hz∼3000Hz,NBFM 调频信号
s
(
t
)
s(t)
s(t) 双边幅度谱有两根离散谱线(
±
20000
H
z
{\pm}20000Hz
±20000Hz)及两个谱峰(范围约为
−
23000
H
z
∼
−
17000
H
z
-23000Hz{\sim}-17000Hz
−23000Hz∼−17000Hz、
17000
H
z
∼
23000
H
z
17000Hz{\sim}23000Hz
17000Hz∼23000Hz)。NBFM 调频信号
s
(
t
)
s(t)
s(t) 的带宽满足公式
B
N
B
F
M
≈
2
f
H
=
6000
H
z
B_{NBFM}\approx2f_H=6000Hz
BNBFM≈2fH=6000Hz,代码详见附录 main_modFM_example2.m
与 mod_fm.m
,设置 beta = 1e-6
即可。
1.5 宽带 FM 信号示例
当
m
(
t
)
m(t)
m(t) 是确知信号时,设
f
m
=
2500
H
z
f_m=2500Hz
fm=2500Hz,
β
=
4
\beta=4
β=4,其他参数与 1.4 节中的确知信号相同,FM 调制效果如下图所示,带宽约为
B
F
M
=
2
(
β
+
1
)
f
m
=
25000
H
z
B_{FM}=2(\beta+1)f_m=25000Hz
BFM=2(β+1)fm=25000Hz,代码详见附录 main_modFM_example1.m
与 mod_fm.m
,设置 fm = 2500, beta = 4
即可。
.当
m
(
t
)
m(t)
m(t) 是确知信号时,设
f
m
=
1000
H
z
f_m=1000Hz
fm=1000Hz,
β
=
15
\beta=15
β=15,其他参数与 1.4 节中的确知信号相同,FM 调制效果如下图所示,带宽约为
B
W
B
F
M
≈
2
Δ
f
=
2
β
f
m
=
30000
H
z
B_{WBFM}\approx2{\Delta}f=2{\beta}f_m=30000Hz
BWBFM≈2Δf=2βfm=30000Hz,代码详见附录 main_modFM_example1.m
与 mod_fm.m
,设置 fm = 1000, beta = 15
即可。
2. FM 解调算法
解调是调制的逆过程,其作用是从接收的已调信号中恢复原基带信号(即调制信号)。FM 解调的方法也分为相干解调和非相干解调,普通的相干解调仅适用于 NBFM 信号,正交相干解调与非相干解调对 NBFM 信号和 WBFM 信号均适用。下面分别用几种不同方法对 1.4 与 1.5 节中确知信号的 FM 调制结果进行解调。
2.1 非相干解调(鉴频器)
微分器和包络检波器是一种最常见的鉴频器。其中微分器的作用是把幅度恒定的调频波
s
F
M
(
t
)
s_{FM}(t)
sFM(t) 变成幅度和频率都随调制信号
m
(
t
)
m(t)
m(t) 变化的调幅调频波
s
d
(
t
)
s_d(t)
sd(t),即:
s
d
(
t
)
=
d
s
F
M
(
t
)
d
t
=
d
{
A
c
o
s
[
ω
c
t
+
K
f
∫
0
t
m
(
τ
)
d
τ
]
}
d
t
=
−
A
[
ω
c
+
K
f
m
(
t
)
]
s
i
n
[
ω
c
t
+
K
f
∫
0
t
m
(
τ
)
d
τ
]
(17)
\begin{aligned} s_d(t)&=\frac{ds_{FM}(t)}{dt}\\[1em] &=\frac{d\left\{ Acos{\left[\omega_ct+{K_f}\int_0^{t}{m(\tau)}d{\tau}\right]}\right\}}{dt}\\[1em] &=-A\left[\omega_c+{K_f}m(t)\right]sin{\left[\omega_ct+{K_f}\int_0^{t}{m(\tau)}d{\tau}\right]}\\[1em] \end{aligned} \tag{17}
sd(t)=dtdsFM(t)=dtd{Acos[ωct+Kf∫0tm(τ)dτ]}=−A[ωc+Kfm(t)]sin[ωct+Kf∫0tm(τ)dτ](17)
包络检波器则将其幅度变化检出并滤去直流,即得解调输出
m
o
(
t
)
=
K
d
K
f
m
(
t
)
m_o(t)={K_d}{K_f}m(t)
mo(t)=KdKfm(t),式中
K
d
K_d
Kd 为鉴频器灵敏度(单位
V
/
(
r
a
d
/
s
)
V/(rad/s)
V/(rad/s))。FM 非相干解调(鉴频器)一般有以下四个步骤,其中第一步是微分器,第二至第四步是包络检波器,与 AM 包络检波的流程一样。
- 第一步:求微分,得到 s d ( t ) s_d(t) sd(t)。
- 第二步:全波整流(对 s ( t ) s(t) s(t) 取绝对值)或半波整流(将 s ( t ) s(t) s(t) 小于 0 0 0 的地方置零)。
- 第三步:低通滤波器滤除高频载波,滤除 2 ω c 2{\omega}_c 2ωc 或 ω c {\omega}_c ωc。
- 第四步:去除直流分量(减去自身均值)。
对 1.4 节中的 FM 信号,设定 f m = 2500 H z f_m=2500Hz fm=2500Hz, β = 1 0 − 6 {\beta}=10^{-6} β=10−6,信噪比 S N R = 150 d B SNR=150dB SNR=150dB,非相干解调效果如下,最大频偏 Δ f = β f m = 0.0025 H z {\Delta}f={\beta}f_m=0.0025Hz Δf=βfm=0.0025Hz 过小,解调效果很差。
对 1.5 节中的 FM 信号,设定 f m = 2500 H z f_m=2500Hz fm=2500Hz, β = 4 {\beta}=4 β=4,信噪比 S N R = 50 d B SNR=50dB SNR=50dB,非相干解调效果如下,解调后幅度放大系数 k = ∣ m ( t ) ∣ ‾ / ∣ m ^ ( t ) ∣ ‾ ≈ 7.61 k=\overline{{\lvert}m(t){\rvert}}/\overline{{\lvert}\hat{m}(t){\rvert}}\approx7.61 k=∣m(t)∣/∣m^(t)∣≈7.61,使用这个系数放大解调信号幅值,然后计算误差,有: ∑ ∣ m ( t i ) − k m ^ ( t i ) ∣ 2 / ∑ ∣ m ( t i ) ∣ 2 ≈ 0.0920 \sqrt{\sum{{\lvert}m(t_i)-k\hat{m}(t_i){\rvert}^2}}/\sqrt{\sum{{\lvert}m(t_i){\rvert}^2}}\approx0.0920 ∑∣m(ti)−km^(ti)∣2/∑∣m(ti)∣2≈0.0920。
对 1.5 节中的 FM 信号,设定 f m = 1000 H z f_m=1000Hz fm=1000Hz, β = 15 {\beta}=15 β=15,信噪比 S N R = 50 d B SNR=50dB SNR=50dB,非相干解调效果如下,解调后幅度放大系数 k = ∣ m ( t ) ∣ ‾ / ∣ m ^ ( t ) ∣ ‾ ≈ 5.07 k=\overline{{\lvert}m(t){\rvert}}/\overline{{\lvert}\hat{m}(t){\rvert}}\approx5.07 k=∣m(t)∣/∣m^(t)∣≈5.07,使用这个系数放大解调信号幅值,然后计算误差,有: ∑ ∣ m ( t i ) − k m ^ ( t i ) ∣ 2 / ∑ ∣ m ( t i ) ∣ 2 ≈ 0.0634 \sqrt{\sum{{\lvert}m(t_i)-k\hat{m}(t_i){\rvert}^2}}/\sqrt{\sum{{\lvert}m(t_i){\rvert}^2}}\approx0.0634 ∑∣m(ti)−km^(ti)∣2/∑∣m(ti)∣2≈0.0634。
非相干解调(鉴频器)的代码详见 lpf_filter.m
、demod_fm_method1.m
和 main_demodFM_example1.m
。
2.2 非相干解调(鉴频器 - 希尔伯特变换法)
鉴频器中的包络检波器可以通过希尔伯特变换法来实现,解调时无需任何载频信息,这样,FM 非相干解调(鉴频器)可以分为以下三个步骤。
- 第一步:求微分,得到 s d ( t ) s_d(t) sd(t)。
- 第二步:计算 s d ( t ) s_d(t) sd(t) 的希尔伯特变换,得到一个复信号(实部为 s d ( t ) s_d(t) sd(t),虚部为其希尔伯特变换结果),对所得复信号取模,即为 s d ( t ) s_d(t) sd(t) 的包络。
- 第三步:去除直流分量(减去自身均值)。
对 1.4 节中的 FM 信号,设定 f m = 2500 H z f_m=2500Hz fm=2500Hz, β = 1 0 − 6 {\beta}=10^{-6} β=10−6,信噪比 S N R = 150 d B SNR=150dB SNR=150dB,非相干解调效果如下,最大频偏 Δ f = β f m = 0.0025 H z {\Delta}f={\beta}f_m=0.0025Hz Δf=βfm=0.0025Hz 过小,解调效果很差。
对 1.5 节中的 FM 信号,设定 f m = 2500 H z f_m=2500Hz fm=2500Hz, β = 4 {\beta}=4 β=4,信噪比 S N R = 50 d B SNR=50dB SNR=50dB,非相干解调效果如下,解调后幅度放大系数 k = ∣ m ( t ) ∣ ‾ / ∣ m ^ ( t ) ∣ ‾ ≈ 4.87 k=\overline{{\lvert}m(t){\rvert}}/\overline{{\lvert}\hat{m}(t){\rvert}}\approx4.87 k=∣m(t)∣/∣m^(t)∣≈4.87,使用这个系数放大解调信号幅值,然后计算误差,有: ∑ ∣ m ( t i ) − k m ^ ( t i ) ∣ 2 / ∑ ∣ m ( t i ) ∣ 2 ≈ 0.0492 \sqrt{\sum{{\lvert}m(t_i)-k\hat{m}(t_i){\rvert}^2}}/\sqrt{\sum{{\lvert}m(t_i){\rvert}^2}}\approx0.0492 ∑∣m(ti)−km^(ti)∣2/∑∣m(ti)∣2≈0.0492。
对 1.5 节中的 FM 信号,设定 f m = 1000 H z f_m=1000Hz fm=1000Hz, β = 15 {\beta}=15 β=15,信噪比 S N R = 50 d B SNR=50dB SNR=50dB,非相干解调效果如下,解调后幅度放大系数 k = ∣ m ( t ) ∣ ‾ / ∣ m ^ ( t ) ∣ ‾ ≈ 3.26 k=\overline{{\lvert}m(t){\rvert}}/\overline{{\lvert}\hat{m}(t){\rvert}}\approx3.26 k=∣m(t)∣/∣m^(t)∣≈3.26,使用这个系数放大解调信号幅值,然后计算误差,有: ∑ ∣ m ( t i ) − k m ^ ( t i ) ∣ 2 / ∑ ∣ m ( t i ) ∣ 2 ≈ 0.0433 \sqrt{\sum{{\lvert}m(t_i)-k\hat{m}(t_i){\rvert}^2}}/\sqrt{\sum{{\lvert}m(t_i){\rvert}^2}}\approx0.0433 ∑∣m(ti)−km^(ti)∣2/∑∣m(ti)∣2≈0.0433。
非相干解调(鉴频器 - 希尔伯特变换法)的代码详见 demod_fm_method2.m
和 main_demodFM_example2.m
。
2.3 相干解调
相干解调时,为了无失真地恢复原基带信号,接收端必须提供一个与调制载波**严格同步(同频同相)**的本地载波(称为相干载波,可使用锁相环技术得到)。普通的相干解调仅适用于 NBFM 信号,其解调框图如下:
其中:
s
N
B
F
M
(
t
)
=
A
c
o
s
(
ω
c
t
)
−
[
A
K
f
∫
0
t
m
(
τ
)
d
τ
]
s
i
n
(
ω
c
t
)
(18)
s_{NBFM}(t)={Acos}{(\omega_ct)}-\left[A{K_f}\int_0^{t}{m(\tau)}d{\tau}\right]sin{(\omega_ct)} \tag{18}
sNBFM(t)=Acos(ωct)−[AKf∫0tm(τ)dτ]sin(ωct)(18)
相干载波
c
(
t
)
=
−
s
i
n
(
ω
c
t
)
c(t)=-sin({\omega_c}t)
c(t)=−sin(ωct),则相乘器的输出为:
s
p
(
t
)
=
−
A
2
s
i
n
(
2
ω
c
t
)
+
[
A
K
f
2
∫
0
t
m
(
τ
)
d
τ
]
⋅
[
1
−
c
o
s
(
2
ω
c
t
)
]
(19)
s_p(t)=-\frac{A}{2}sin{(2{\omega_c}t)}+\left[\frac{AK_f}{2}\int_0^{t}{m(\tau)}d{\tau}\right]\cdot\left[{1-cos(2{\omega_c}t)}\right] \tag{19}
sp(t)=−2Asin(2ωct)+[2AKf∫0tm(τ)dτ]⋅[1−cos(2ωct)](19)
经低通滤波器取出其低频分量:
s
d
(
t
)
=
A
K
f
2
∫
0
t
m
(
τ
)
d
τ
(20)
s_d(t)=\frac{AK_f}{2}\int_0^{t}{m(\tau)}d{\tau} \tag{20}
sd(t)=2AKf∫0tm(τ)dτ(20)
再经微分器,即得解调输出:
m
o
(
t
)
=
A
K
f
2
m
(
t
)
(21)
m_o(t)=\frac{AK_f}{2}m(t) \tag{21}
mo(t)=2AKfm(t)(21)
可见,相干解调可以恢复原调制信号,这种解调方法与线性调制中的相干解调一样,要求本地载波与调制载波同步,否则将使解调信号失真。NBFM 相干解调可以分为以下几步:
- 第一步:乘以相干载波(即乘以 − s i n ( ω c t + ϕ 0 ) -sin({\omega_c}t+\phi_0) −sin(ωct+ϕ0)),注意载波初始相位。
- 第二步:低通滤波滤除高频载波,滤除 2 ω c 2\omega_c 2ωc。
- 第三步:求微分。
对 1.4 节中的 FM 信号,设定 f m = 2500 H z f_m=2500Hz fm=2500Hz, β = 1 0 − 6 {\beta}=10^{-6} β=10−6,信噪比 S N R = 150 d B SNR=150dB SNR=150dB,相干解调效果如下,解调后幅度放大系数 k = ∣ m ( t ) ∣ ‾ / ∣ m ^ ( t ) ∣ ‾ ≈ 222.30 k=\overline{{\lvert}m(t){\rvert}}/\overline{{\lvert}\hat{m}(t){\rvert}}\approx222.30 k=∣m(t)∣/∣m^(t)∣≈222.30,使用这个系数放大解调信号幅值,然后计算误差,有: ∑ ∣ m ( t i ) − k m ^ ( t i ) ∣ 2 / ∑ ∣ m ( t i ) ∣ 2 ≈ 0.1313 \sqrt{\sum{{\lvert}m(t_i)-k\hat{m}(t_i){\rvert}^2}}/\sqrt{\sum{{\lvert}m(t_i){\rvert}^2}}\approx0.1313 ∑∣m(ti)−km^(ti)∣2/∑∣m(ti)∣2≈0.1313。更改相干载波的初始相位为 ϕ 0 = π / 4 , π / 2 {\phi_0}=\pi/4,\pi/2 ϕ0=π/4,π/2,或者更改相干载波的中心频率为 0.8 f c , 1.2 f c 0.8f_c,1.2f_c 0.8fc,1.2fc 后,解调效果变差,说明这种方法对相干载波同频同相的要求较高,鲁棒性不够强悍,可使用锁相环技术来改善这一缺点。
对 1.5 节中的 FM 信号,设定 f m = 2500 H z f_m=2500Hz fm=2500Hz, β = 4 {\beta}=4 β=4,信噪比 S N R = 50 d B SNR=50dB SNR=50dB,相干解调效果如下,可知这种方法不适用于非窄带调频信号。
对 1.5 节中的 FM 信号,设定 f m = 1000 H z f_m=1000Hz fm=1000Hz, β = 15 {\beta}=15 β=15,信噪比 S N R = 50 d B SNR=50dB SNR=50dB,相干解调效果如下,可知这种方法不适用于非窄带调频信号。
相干解调的代码详见 lpf_filter.m
、demod_fm_method3.m
和 main_demodFM_example3.m
。
2.4 数字正交解调
数字正交解调也属于相干解调的一种,但这种方法具有较强的抗载频失配能力,不要求相干载波严格的同频同相。FM 数字正交解调一般有以下四个步骤:
- 第一步:乘以正交相干载波得到 s I ( t ) {s_I}(t) sI(t) 与 s Q ( t ) {s_Q}(t) sQ(t),即 s I ( t ) = s ( t ) c o s ( ω c t + ϕ 0 ) {s_I}(t)=s(t)cos({\omega_ct}+{\phi_0}) sI(t)=s(t)cos(ωct+ϕ0), s Q ( t ) = − s ( t ) s i n ( ω c t + ϕ 0 ) {s_Q}(t)=-s(t)sin({\omega_ct}+{\phi_0}) sQ(t)=−s(t)sin(ωct+ϕ0)。
- 第二步:低通滤波器滤除 s I ( t ) {s_I}(t) sI(t) 与 s Q ( t ) {s_Q}(t) sQ(t) 中的高频分量。
- 第三步:通过反正切计算相位 Φ ( t ) = a t a n [ s Q ( t ) s I ( t ) ] = k ∫ 0 t m ( τ ) d τ \Phi(t)=atan\left[\frac{s_Q(t)}{s_I(t)}\right]=k\int_0^{t}{m(\tau)}d{\tau} Φ(t)=atan[sI(t)sQ(t)]=k∫0tm(τ)dτ。
- 第四步:对相位进行差分,得到解调结果 m o ( t ) m_o(t) mo(t)。
反正切运算在硬件中难以实现,通过一阶近似,上面第三步与第四步可合并为以下式子:
m
o
(
t
)
=
s
I
(
n
−
1
)
s
Q
(
n
)
−
s
I
(
n
)
s
Q
(
n
−
1
)
s
I
2
(
n
)
+
s
Q
2
(
n
)
(22)
m_o(t)=\frac{{s_I(n-1)}{s_Q(n)}-{s_I(n)}{s_Q(n-1)}}{{s_I^2(n)}+{s_Q^2(n)}} \tag{22}
mo(t)=sI2(n)+sQ2(n)sI(n−1)sQ(n)−sI(n)sQ(n−1)(22)
对 1.4 节中的 FM 信号,设定
f
m
=
2500
H
z
f_m=2500Hz
fm=2500Hz,
β
=
1
0
−
6
{\beta}=10^{-6}
β=10−6,信噪比
S
N
R
=
150
d
B
SNR=150dB
SNR=150dB,数字正交解调效果如下,解调后幅度放大系数
k
=
∣
m
(
t
)
∣
‾
/
∣
m
^
(
t
)
∣
‾
≈
17782852.62
k=\overline{{\lvert}m(t){\rvert}}/\overline{{\lvert}\hat{m}(t){\rvert}}\approx17782852.62
k=∣m(t)∣/∣m^(t)∣≈17782852.62,使用这个系数放大解调信号幅值,然后计算误差,有:
∑
∣
m
(
t
i
)
−
k
m
^
(
t
i
)
∣
2
/
∑
∣
m
(
t
i
)
∣
2
≈
0.1318
\sqrt{\sum{{\lvert}m(t_i)-k\hat{m}(t_i){\rvert}^2}}/\sqrt{\sum{{\lvert}m(t_i){\rvert}^2}}\approx0.1318
∑∣m(ti)−km^(ti)∣2/∑∣m(ti)∣2≈0.1318。
对 1.5 节中的 FM 信号,设定 f m = 2500 H z f_m=2500Hz fm=2500Hz, β = 4 {\beta}=4 β=4,信噪比 S N R = 50 d B SNR=50dB SNR=50dB,数字正交解调效果如下,解调后幅度放大系数 k = ∣ m ( t ) ∣ ‾ / ∣ m ^ ( t ) ∣ ‾ ≈ 4.48 k=\overline{{\lvert}m(t){\rvert}}/\overline{{\lvert}\hat{m}(t){\rvert}}\approx4.48 k=∣m(t)∣/∣m^(t)∣≈4.48,使用这个系数放大解调信号幅值,然后计算误差,有: ∑ ∣ m ( t i ) − k m ^ ( t i ) ∣ 2 / ∑ ∣ m ( t i ) ∣ 2 ≈ 0.0391 \sqrt{\sum{{\lvert}m(t_i)-k\hat{m}(t_i){\rvert}^2}}/\sqrt{\sum{{\lvert}m(t_i){\rvert}^2}}\approx0.0391 ∑∣m(ti)−km^(ti)∣2/∑∣m(ti)∣2≈0.0391。
对 1.5 节中的 FM 信号,设定 f m = 1000 H z f_m=1000Hz fm=1000Hz, β = 15 {\beta}=15 β=15,信噪比 S N R = 50 d B SNR=50dB SNR=50dB,数字正交解调效果如下,解调后幅度放大系数 k = ∣ m ( t ) ∣ ‾ / ∣ m ^ ( t ) ∣ ‾ ≈ 2.99 k=\overline{{\lvert}m(t){\rvert}}/\overline{{\lvert}\hat{m}(t){\rvert}}\approx2.99 k=∣m(t)∣/∣m^(t)∣≈2.99,使用这个系数放大解调信号幅值,然后计算误差,有: ∑ ∣ m ( t i ) − k m ^ ( t i ) ∣ 2 / ∑ ∣ m ( t i ) ∣ 2 ≈ 0.0163 \sqrt{\sum{{\lvert}m(t_i)-k\hat{m}(t_i){\rvert}^2}}/\sqrt{\sum{{\lvert}m(t_i){\rvert}^2}}\approx0.0163 ∑∣m(ti)−km^(ti)∣2/∑∣m(ti)∣2≈0.0163。
数字正交解调的代码详见 lpf_filter.m
、demod_fm_method4.m
和 main_demodFM_example4.m
。更改相干载波的初始相位为
ϕ
0
=
π
/
4
,
π
/
2
{\phi_0}=\pi/4,\pi/2
ϕ0=π/4,π/2,或者更改相干载波的中心频率为
0.8
f
c
,
1.2
f
c
0.8f_c,1.2f_c
0.8fc,1.2fc 后,解调效果依然比较好,说明这种方法具有较好的抗载频失配能力。
3. FM 仿真(MATLAB Communications Toolbox)
MATLAB 的 Communications Toolbox 中提供了 FM 调制函数 fmmod,高斯白噪声函数 awgn,以及 FM 解调函数 fmdemod,可以很方便地完成 FM 信号仿真。使用这三个函数实现上面 1.4 节中确知信号 m ( t ) m(t) m(t) 的 FM 调制解调,调制后加噪声的效果如下:
解调效果如下:
解调信号与调制信号波形基本重回,计算误差,有:
∑
∣
m
(
t
i
)
−
m
^
(
t
i
)
∣
2
/
∑
∣
m
(
t
i
)
∣
2
≈
0.2397
\sqrt{\sum{{\lvert}m(t_i)-\hat{m}(t_i){\rvert}^2}}/\sqrt{\sum{{\lvert}m(t_i){\rvert}^2}}\approx0.2397
∑∣m(ti)−m^(ti)∣2/∑∣m(ti)∣2≈0.2397。代码详见附录 main_CommFM_example.m
。
参考资料
[1] 楼才义,徐建良,杨小牛.软件无线电原理与应用[M].电子工业出版社,2014.
[2] 樊昌信,曹丽娜.通信原理.第7版[M].国防工业出版社,2012.
[3] 刘学勇.详解MATLAB/Simulink通信系统建模与仿真[M].电子工业出版社,2011.
[4] 王丽娜,王兵.卫星通信系统.第2版[M].国防工业出版社,2014.
[5] 博客园 - 两种频率调制(FM)方法的MATLAB实现。
[6] CSDN - 数字信号处理基础----FM的调制与解调。
附录代码
附.1 文件 mod_fm.m
function [ sig_fm, deltaf ] = mod_fm(fc, beta, fs, mt, t, W, A)
% MOD_FM FM 调频
% 输入参数:
% fc 载波中心频率
% beta 调频指数/调制指数
% fs 信号采样率
% mt 调制信号
% t 采样时间
% W 基带信号带宽(最高频率)
% A 载波恒定振幅
% 输出参数:
% sig_fm 调频(FM)实信号
% deltaf 最大频偏
% @author 木三百川
% 计算调频灵敏度及最大频偏
kf = beta*W/max(abs(mt));
deltaf = beta*W;
% 计算调制信号积分
int_mt = cumtrapz(t,mt);
% 生成信号
sig_fm = A*cos(2*pi*fc*t+2*pi*kf*int_mt); % FM调频信号
% 绘图
nfft = length(sig_fm);
freq = (-nfft/2:nfft/2-1).'*(fs/nfft);
figure;set(gcf,'color','w');
plot_length = min(500, length(sig_fm));
subplot(3,2,1);
plot(t(1:plot_length), mt(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('调制信号m(t)');
subplot(3,2,2);
plot(freq, 10*log10(fftshift(abs(fft(mt,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('频率/hz');ylabel('幅度/dB');title('调制信号m(t)双边幅度谱');
subplot(3,2,3);
plot(t(1:plot_length), int_mt(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('调制信号m(t)积分');
subplot(3,2,4);
plot(freq, 10*log10(fftshift(abs(fft(int_mt,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('频率/hz');ylabel('幅度/dB');title('调制信号m(t)积分双边幅度谱');
subplot(3,2,5);
plot(t(1:plot_length), sig_fm(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('FM调频信号s(t)');
subplot(3,2,6);
plot(freq, 10*log10(fftshift(abs(fft(sig_fm,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('频率/hz');ylabel('幅度/dB');title('FM调频信号s(t)双边幅度谱');
end
附.2 文件 main_modFM_example1.m
clc;
clear;
close all;
% FM 调制仿真(调制信号为确知信号)
% @author 木三百川
% 调制参数
A = 1; % 载波恒定振幅
fm = 2500; % 调制信号参数
beta = 1e-6; % 调频指数/调制指数
fc = 20000; % 载波频率
fs = 8*fc; % 采样率
total_time = 2; % 仿真时长,单位:秒
% 采样时间
t = 0:1/fs:total_time-1/fs;
% 调制信号为确知信号
mt = sin(2*pi*fm*t)+cos(pi*fm*t);
W = fm;
% FM 调制
[ sig_fm, deltaf ] = mod_fm(fc, beta, fs, mt, t, W, A);
fprintf('最大频偏 deltaf = %.6f Hz.\n', deltaf);
附.3 文件 main_modFM_example2.m
clc;
clear;
close all;
% FM 调制仿真(调制信号为随机信号)
% @author 木三百川
% 调制参数
A = 1; % 载波恒定振幅
fH = 3000; % 基带信号带宽
beta = 1e-6; % 调频指数/调制指数
fc = 20000; % 载波频率
fs = 8*fc; % 采样率
total_time = 2; % 仿真时长,单位:秒
% 采样时间
t = 0:1/fs:total_time-1/fs;
% 调制信号为随机信号
mt = randn(size(t));
b = fir1(512, fH/(fs/2), 'low');
mt = filter(b,1,mt);
mt = mt - mean(mt);
W = fH;
% FM 调制
[ sig_fm, deltaf ] = mod_fm(fc, beta, fs, mt, t, W, A);
fprintf('最大频偏 deltaf = %.6f Hz.\n', deltaf);
附.4 文件 lpf_filter.m
function sig_lpf = lpf_filter(sig_data, cutfre)
% LPF_FILTER 自定义理想低通滤波器
% 输入参数:
% sig_data 待滤波数据
% cutfre 截止频率,范围 (0,1)
% 输出参数:
% sig_lpf 低通滤波结果
% @author 木三百川
nfft = length(sig_data);
lidx = round(nfft/2-cutfre*nfft/2);
ridx = nfft - lidx;
sig_fft_lpf = fftshift(fft(sig_data));
sig_fft_lpf([1:lidx,ridx:nfft]) = 0;
sig_lpf = real(ifft(fftshift(sig_fft_lpf)));
end
附.5 文件 demod_fm_method1.m
function [ sig_fm_demod ] = demod_fm_method1(sig_fm_receive, fc, fs, t)
% DEMOD_FM_METHOD1 FM 非相干解调(鉴频器)
% 输入参数:
% sig_fm_receive FM 接收信号,行向量
% fc 载波中心频率
% fs 信号采样率
% t 采样时间
% 输出参数:
% sig_fm_demod 解调结果,与 sig_fm_receive 等长
% @author 木三百川
% 第一步:求微分
sig_dfm = [diff(sig_fm_receive),0];
% 第二步:全波整流
sig_fm_abs = abs(sig_dfm);
% 第三步:低通滤波
sig_fm_lpf = lpf_filter(sig_fm_abs, fc/2/(fs/2));
% 第四步:去除直流分量
sig_fm_demod = sig_fm_lpf - mean(sig_fm_lpf);
% 绘图
nfft = length(sig_fm_abs);
freq = (-nfft/2:nfft/2-1).'*(fs/nfft);
figure;set(gcf,'color','w');
plot_length = min(500, length(sig_fm_abs));
subplot(3,2,1);
plot(t(1:plot_length), sig_fm_abs(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('全波整流结果');
subplot(3,2,2);
plot(freq, 10*log10(fftshift(abs(fft(sig_fm_abs,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('频率/hz');ylabel('幅度/dB');title('全波整流结果双边幅度谱');
subplot(3,2,3);
plot(t(1:plot_length), sig_fm_lpf(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('低通滤波结果');
subplot(3,2,4);
plot(freq, 10*log10(fftshift(abs(fft(sig_fm_lpf,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('频率/hz');ylabel('幅度/dB');title('低通滤波结果双边幅度谱');
subplot(3,2,5);
plot(t(1:plot_length), sig_fm_demod(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('(去除直流)解调结果');
subplot(3,2,6);
plot(freq, 10*log10(fftshift(abs(fft(sig_fm_demod,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('频率/hz');ylabel('幅度/dB');title('(去除直流)解调结果双边幅度谱');
end
附.6 文件 main_demodFM_example1.m
clc;
clear;
close all;
% FM 解调仿真(调制信号为确知信号,非相干解调)
% @author 木三百川
% 调制参数
A = 1; % 载波恒定振幅
fm = 1000; % 调制信号参数
beta = 15; % 调频指数/调制指数
fc = 20000; % 载波频率
fs = 8*fc; % 采样率
total_time = 2; % 仿真时长,单位:秒
% 采样时间
t = 0:1/fs:total_time-1/fs;
% 调制信号为确知信号
mt = sin(2*pi*fm*t)+cos(pi*fm*t);
W = fm;
% FM 调制
[ sig_fm_send, deltaf ] = mod_fm(fc, beta, fs, mt, t, W, A);
fprintf('最大频偏 deltaf = %.6f Hz.\n', deltaf);
% 加噪声
snr = 50; % 信噪比
sig_fm_receive = awgn(sig_fm_send, snr, 'measured');
% 非相干解调
[ sig_fm_demod ] = demod_fm_method1(sig_fm_receive, fc, fs, t);
% 绘图
nfft = length(sig_fm_receive);
freq = (-nfft/2:nfft/2-1).'*(fs/nfft);
figure;set(gcf,'color','w');
plot_length = min(500, length(sig_fm_receive));
subplot(1,2,1);
plot(t(1:plot_length), sig_fm_receive(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('FM接收信号');
subplot(1,2,2);
plot(freq, 10*log10(fftshift(abs(fft(sig_fm_receive,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('频率/hz');ylabel('幅度/dB');title('FM接收信号双边幅度谱');
coef = mean(abs(mt))/mean(abs(sig_fm_demod));
fprintf('norm(调制信号 - %.2f * 解调信号)/norm(调制信号) = %.4f.\n', coef, norm(mt-coef*sig_fm_demod)/norm(mt));
figure;set(gcf,'color','w');
plot(t(1:plot_length), mt(1:plot_length));xlim([t(1),t(plot_length)]);
hold on;
plot(t(1:plot_length), coef*sig_fm_demod(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('解调效果');
legend('调制信号','解调信号(放大后)');
附.7 文件 demod_fm_method2.m
function [ sig_fm_demod ] = demod_fm_method2(sig_fm_receive, fs, t)
% DEMOD_FM_METHOD2 FM 非相干解调(鉴频器 - 希尔伯特变换法)
% 输入参数:
% sig_fm_receive FM 接收信号,行向量
% fs 信号采样率
% t 采样时间
% 输出参数:
% sig_fm_demod 解调结果,与 sig_fm_receive 等长
% @author 木三百川
% 第一步:求微分
sig_dfm = [diff(sig_fm_receive),0];
% 第二步:计算信号包络
sig_fm_envelope = abs(hilbert(sig_dfm));
% 第三步:去除直流分量
sig_fm_demod = sig_fm_envelope - mean(sig_fm_envelope);
% 绘图
nfft = length(sig_fm_receive);
freq = (-nfft/2:nfft/2-1).'*(fs/nfft);
figure;set(gcf,'color','w');
plot_length = min(500, length(sig_fm_receive));
subplot(3,2,1);
plot(t(1:plot_length), sig_dfm(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('求微分结果');
subplot(3,2,2);
plot(freq, 10*log10(fftshift(abs(fft(sig_dfm,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('频率/hz');ylabel('幅度/dB');title('求微分结果双边幅度谱');
subplot(3,2,3);
plot(t(1:plot_length), sig_fm_envelope(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('计算信号包络结果');
subplot(3,2,4);
plot(freq, 10*log10(fftshift(abs(fft(sig_fm_envelope,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('频率/hz');ylabel('幅度/dB');title('计算信号包络结果双边幅度谱');
subplot(3,2,5);
plot(t(1:plot_length), sig_fm_demod(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('(去除直流)解调结果');
subplot(3,2,6);
plot(freq, 10*log10(fftshift(abs(fft(sig_fm_demod,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('频率/hz');ylabel('幅度/dB');title('(去除直流)解调结果双边幅度谱');
end
附.8 文件 main_demodFM_example2.m
clc;
clear;
close all;
% FM 解调仿真(调制信号为确知信号,非相干解调/希尔伯特变换法)
% @author 木三百川
% 调制参数
A = 1; % 载波恒定振幅
fm = 1000; % 调制信号参数
beta = 15; % 调频指数/调制指数
fc = 20000; % 载波频率
fs = 8*fc; % 采样率
total_time = 2; % 仿真时长,单位:秒
% 采样时间
t = 0:1/fs:total_time-1/fs;
% 调制信号为确知信号
mt = sin(2*pi*fm*t)+cos(pi*fm*t);
W = fm;
% FM 调制
[ sig_fm_send, deltaf ] = mod_fm(fc, beta, fs, mt, t, W, A);
fprintf('最大频偏 deltaf = %.6f Hz.\n', deltaf);
% 加噪声
snr = 50; % 信噪比
sig_fm_receive = awgn(sig_fm_send, snr, 'measured');
% 非相干解调
[ sig_fm_demod ] = demod_fm_method2(sig_fm_receive, fs, t);
% 绘图
nfft = length(sig_fm_receive);
freq = (-nfft/2:nfft/2-1).'*(fs/nfft);
figure;set(gcf,'color','w');
plot_length = min(500, length(sig_fm_receive));
subplot(1,2,1);
plot(t(1:plot_length), sig_fm_receive(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('FM接收信号');
subplot(1,2,2);
plot(freq, 10*log10(fftshift(abs(fft(sig_fm_receive,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('频率/hz');ylabel('幅度/dB');title('FM接收信号双边幅度谱');
coef = mean(abs(mt))/mean(abs(sig_fm_demod));
fprintf('norm(调制信号 - %.2f * 解调信号)/norm(调制信号) = %.4f.\n', coef, norm(mt-coef*sig_fm_demod)/norm(mt));
figure;set(gcf,'color','w');
plot(t(1:plot_length), mt(1:plot_length));xlim([t(1),t(plot_length)]);
hold on;
plot(t(1:plot_length), coef*sig_fm_demod(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('解调效果');
legend('调制信号','解调信号(放大后)');
附.9 文件 demod_fm_method3.m
function [ sig_fm_demod ] = demod_fm_method3(sig_fm_receive, fc, fs, t, phi0)
% DEMOD_FM_METHOD3 FM 相干解调
% 输入参数:
% sig_fm_receive FM 接收信号,行向量
% fc 载波中心频率
% fs 信号采样率
% t 采样时间
% phi0 相干载波初始相位
% 输出参数:
% sig_fm_demod 解调结果,与 sig_fm_receive 等长
% @author 木三百川
% 第一步:乘以相干载波
sig_fm_ct = -sig_fm_receive.*sin(2*pi*fc*t+phi0);
% 第二步:低通滤波
sig_fm_lpf = lpf_filter(sig_fm_ct, fc/(fs/2));
% 第三步:求微分
sig_fm_demod = [diff(sig_fm_lpf),0]*fs;
% 绘图
nfft = length(sig_fm_receive);
freq = (-nfft/2:nfft/2-1).'*(fs/nfft);
figure;set(gcf,'color','w');
plot_length = min(500, length(sig_fm_receive));
subplot(3,2,1);
plot(t(1:plot_length), sig_fm_ct(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('乘以相干载波结果');
subplot(3,2,2);
plot(freq, 10*log10(fftshift(abs(fft(sig_fm_ct,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('频率/hz');ylabel('幅度/dB');title('乘以相干载波结果双边幅度谱');
subplot(3,2,3);
plot(t(1:plot_length), sig_fm_lpf(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('低通滤波结果');
subplot(3,2,4);
plot(freq, 10*log10(fftshift(abs(fft(sig_fm_lpf,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('频率/hz');ylabel('幅度/dB');title('低通滤波结果双边幅度谱');
subplot(3,2,5);
plot(t(1:plot_length), sig_fm_demod(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('(求微分)解调结果');
subplot(3,2,6);
plot(freq, 10*log10(fftshift(abs(fft(sig_fm_demod,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('频率/hz');ylabel('幅度/dB');title('(求微分)解调结果双边幅度谱');
end
附.10 文件 main_demodFM_example3.m
clc;
clear;
close all;
% FM 解调仿真(调制信号为确知信号,相干解调)
% @author 木三百川
% 调制参数
A = 1; % 载波恒定振幅
fm = 2500; % 调制信号参数
beta = 1e-6; % 调频指数/调制指数
fc = 20000; % 载波频率
fs = 8*fc; % 采样率
total_time = 2; % 仿真时长,单位:秒
% 采样时间
t = 0:1/fs:total_time-1/fs;
% 调制信号为确知信号
mt = sin(2*pi*fm*t)+cos(pi*fm*t);
W = fm;
% FM 调制
[ sig_fm_send, deltaf ] = mod_fm(fc, beta, fs, mt, t, W, A);
fprintf('最大频偏 deltaf = %.6f Hz.\n', deltaf);
% 加噪声
snr = 150; % 信噪比
sig_fm_receive = awgn(sig_fm_send, snr, 'measured');
% 相干解调
ini_phase = 0;
[ sig_fm_demod ] = demod_fm_method3(sig_fm_receive, fc, fs, t, ini_phase);
% 绘图
nfft = length(sig_fm_receive);
freq = (-nfft/2:nfft/2-1).'*(fs/nfft);
figure;set(gcf,'color','w');
plot_length = min(500, length(sig_fm_receive));
subplot(1,2,1);
plot(t(1:plot_length), sig_fm_receive(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('FM接收信号');
subplot(1,2,2);
plot(freq, 10*log10(fftshift(abs(fft(sig_fm_receive,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('频率/hz');ylabel('幅度/dB');title('FM接收信号双边幅度谱');
coef = mean(abs(mt))/mean(abs(sig_fm_demod));
fprintf('norm(调制信号 - %.2f * 解调信号)/norm(调制信号) = %.4f.\n', coef, norm(mt-coef*sig_fm_demod)/norm(mt));
figure;set(gcf,'color','w');
plot(t(1:plot_length), mt(1:plot_length));xlim([t(1),t(plot_length)]);
hold on;
plot(t(1:plot_length), coef*sig_fm_demod(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('解调效果');
legend('调制信号','解调信号(放大后)');
附.11 文件 demod_fm_method4.m
function [ sig_fm_demod ] = demod_fm_method4(sig_fm_receive, fc, fs, t, phi0)
% DEMOD_FM_METHOD4 FM 数字正交解调/相干解调
% 输入参数:
% sig_fm_receive FM 接收信号,行向量
% fc 载波中心频率
% fs 信号采样率
% t 采样时间
% phi0 相干载波初始相位
% 输出参数:
% sig_fm_demod 解调结果,与 sig_fm_receive 等长
% @author 木三百川
% 第一步:乘以正交相干载波
sig_fm_i = sig_fm_receive.*cos(2*pi*fc*t+phi0);
sig_fm_q = -sig_fm_receive.*sin(2*pi*fc*t+phi0);
% 第二步:低通滤波
sig_fm_i_lpf = lpf_filter(sig_fm_i, fc/(fs/2));
sig_fm_q_lpf = lpf_filter(sig_fm_q, fc/(fs/2));
% 第三步:计算相位
Pt = unwrap(atan2(sig_fm_q_lpf, sig_fm_i_lpf));
% 第四步:计算差分
sig_fm_demod = [diff(Pt),0];
% % 合并第三步与第四步:反正切近似
% sig_fm_demod = (sig_fm_i_lpf(1:end-1).*sig_fm_q_lpf(2:end)-sig_fm_i_lpf(2:end).* ...
% sig_fm_q_lpf(1:end-1))./(sig_fm_i_lpf(2:end).^2+sig_fm_q_lpf(2:end).^2);
% sig_fm_demod = [sig_fm_demod, sig_fm_demod(end)];
% 绘图
nfft = length(sig_fm_receive);
freq = (-nfft/2:nfft/2-1).'*(fs/nfft);
figure;set(gcf,'color','w');
plot_length = min(500, length(sig_fm_receive));
subplot(2,2,1);
plot(t(1:plot_length), sig_fm_i(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('乘以正交相干载波 I 路结果');
subplot(2,2,2);
plot(freq, 10*log10(fftshift(abs(fft(sig_fm_i,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('频率/hz');ylabel('幅度/dB');title('乘以正交相干载波 I 路结果双边幅度谱');
subplot(2,2,3);
plot(t(1:plot_length), sig_fm_q(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('乘以正交相干载波 Q 路结果');
subplot(2,2,4);
plot(freq, 10*log10(fftshift(abs(fft(sig_fm_q,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('频率/hz');ylabel('幅度/dB');title('乘以正交相干载波 Q 路结果双边幅度谱');
figure;set(gcf,'color','w');
subplot(2,2,1);
plot(t(1:plot_length), sig_fm_i_lpf(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('低通滤波 I 路结果');
subplot(2,2,2);
plot(freq, 10*log10(fftshift(abs(fft(sig_fm_i_lpf,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('频率/hz');ylabel('幅度/dB');title('低通滤波 I 路结果双边幅度谱');
subplot(2,2,3);
plot(t(1:plot_length), sig_fm_q_lpf(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('低通滤波 Q 路结果');
subplot(2,2,4);
plot(freq, 10*log10(fftshift(abs(fft(sig_fm_q_lpf,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('频率/hz');ylabel('幅度/dB');title('低通滤波 Q 路结果双边幅度谱');
figure;set(gcf,'color','w');
subplot(2,2,1);
plot(t(1:plot_length), Pt(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('计算相位结果');
subplot(2,2,2);
plot(freq, 10*log10(fftshift(abs(fft(Pt,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('频率/hz');ylabel('幅度/dB');title('计算相位结果双边幅度谱');
subplot(2,2,3);
plot(t(1:plot_length), sig_fm_demod(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('(计算差分)解调结果');
subplot(2,2,4);
plot(freq, 10*log10(fftshift(abs(fft(sig_fm_demod,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('频率/hz');ylabel('幅度/dB');title('(计算差分)解调结果双边幅度谱');
end
附.12 文件 main_demodFM_example4.m
clc;
clear;
close all;
% FM 解调仿真(调制信号为确知信号,数字正交解调)
% @author 木三百川
% 调制参数
A = 1; % 载波恒定振幅
fm = 1000; % 调制信号参数
beta = 15; % 调频指数/调制指数
fc = 20000; % 载波频率
fs = 8*fc; % 采样率
total_time = 2; % 仿真时长,单位:秒
% 采样时间
t = 0:1/fs:total_time-1/fs;
% 调制信号为确知信号
mt = sin(2*pi*fm*t)+cos(pi*fm*t);
W = fm;
% FM 调制
[ sig_fm_send, deltaf ] = mod_fm(fc, beta, fs, mt, t, W, A);
fprintf('最大频偏 deltaf = %.6f Hz.\n', deltaf);
% 加噪声
snr = 50; % 信噪比
sig_fm_receive = awgn(sig_fm_send, snr, 'measured');
% 相干解调
ini_phase = 0;
[ sig_fm_demod ] = demod_fm_method4(sig_fm_receive, fc, fs, t, ini_phase);
% 绘图
nfft = length(sig_fm_receive);
freq = (-nfft/2:nfft/2-1).'*(fs/nfft);
figure;set(gcf,'color','w');
plot_length = min(500, length(sig_fm_receive));
subplot(1,2,1);
plot(t(1:plot_length), sig_fm_receive(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('FM接收信号');
subplot(1,2,2);
plot(freq, 10*log10(fftshift(abs(fft(sig_fm_receive,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('频率/hz');ylabel('幅度/dB');title('FM接收信号双边幅度谱');
coef = mean(abs(mt))/mean(abs(sig_fm_demod));
fprintf('norm(调制信号 - %.2f * 解调信号)/norm(调制信号) = %.4f.\n', coef, norm(mt-coef*sig_fm_demod)/norm(mt));
figure;set(gcf,'color','w');
plot(t(1:plot_length), mt(1:plot_length));xlim([t(1),t(plot_length)]);
hold on;
plot(t(1:plot_length), coef*sig_fm_demod(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('解调效果');
legend('调制信号','解调信号(放大后)');
附.13 文件 main_CommFM_example.m
clc;
clear;
close all;
% FM 调制解调仿真(使用Communications Toolbox工具箱)
% @author 木三百川
% 调制参数
A = 1; % 载波恒定振幅
fm = 2500; % 调制信号参数
beta = 1e-6; % 调频指数/调制指数
fc = 20000; % 载波频率
fs = 8*fc; % 采样率
total_time = 2; % 仿真时长,单位:秒
% 采样时间
t = 0:1/fs:total_time-1/fs;
% 调制信号为确知信号
mt = sin(2*pi*fm*t)+cos(pi*fm*t);
% FM 调制
freqdev = beta*fm;
ini_phase = 0;
sig_fm_send = A*fmmod(mt, fc, fs, freqdev, ini_phase);
% 加噪声
snr = 150; % 信噪比
sig_fm_receive = awgn(sig_fm_send, snr, 'measured');
% FM 解调
[ sig_fm_demod ] = fmdemod(sig_fm_receive, fc, fs, freqdev, ini_phase);
% 绘图
nfft = length(sig_fm_receive);
freq = (-nfft/2:nfft/2-1).'*(fs/nfft);
figure;set(gcf,'color','w');
plot_length = min(500, length(sig_fm_receive));
subplot(1,2,1);
plot(t(1:plot_length), sig_fm_receive(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('FM接收信号');
subplot(1,2,2);
plot(freq, 10*log10(fftshift(abs(fft(sig_fm_receive,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('频率/hz');ylabel('幅度/dB');title('FM接收信号双边幅度谱');
figure;set(gcf,'color','w');
plot(t(1:plot_length), mt(1:plot_length));xlim([t(1),t(plot_length)]);
hold on;
plot(t(1:plot_length), sig_fm_demod(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('解调效果');
legend('调制信号','解调信号');
fprintf('norm(调制信号 - 解调信号)/norm(调制信号) = %.4f.\n', norm(mt-sig_fm_demod)/norm(mt));