FM调制解调以及MATLAB实现

89 阅读 0 评论 59 点赞

我是靠谱客的博主健忘篮球，这篇文章主要介绍FM调制解调以及MATLAB实现，现在分享给大家，希望可以做个参考。

首先FM调制时属于非线性调制中的一种，常见的非线性调制时FM(频率调制)，PM(相位)调制
FM调制：就是利用调制信号的变化去改变载波的频率，幅度保持不变，即载波信号的频率随我们调制信号的幅度变化而变化。
首先我们先介绍一些通信原理上的基本得概念：

现在FM调制就是我们将我们的基带信号替换位其φt,但是我们会在其基带信号的前面乘以一个调频系数，也叫调频灵敏度。这就通过我们基带信号的振幅，改变其载波的频率，这样就是实现了频率调制。下面我们假设基带信号为一个余弦波：

上述就是其FM 调制的数学式子，这就实现了其FM 调制

这里说明上图中的积分转化为求和，因为积分是我们的连续域的，而我们的这都是数字信号，都是连续域，这就相当于我们的信号与系统中的傅里叶变化对应连续域，Z变化对应离散域。所以连续域的积分对应到离散域就是求和。

IQ调制：现代通信系统为了使频谱利用率更高，用了许多种矢量调制,IQ调制中，两个正交信号（频率相同，相位相差90 °的载波，一般用Sin和Cos表示）与I（In-Phase，同相分量）、Q（Quadrature Phase，正交分量）两路信号分别调制后一起发射，从而提高频谱利用率。

所以接下来我们要将上式转为IQ两路，然后利用天线进行发送

我们I路信号为：cos ，图中的cos(wt) 为载波信号

我们Q路信号为：sin , 图中的sin(wt) 为载波信号

FM解调：

相干解调：相干解调也叫同步检波，它适用于所有线性调制信号的解调。实现相干解调的关键是接收端要恢复出一个与调制载波严格同步的相干载波。相干解调是指利用乘法器，输入一路与载频相干（同频同相）的参考信号与载频相乘。

非相干解调：通信接收端从已调高频信号中恢复出原始数字基带信号时，采用的非相干解调方式，相对于相干解调方式，是指不需要提取载波信息的一种解调方法。

所谓相干，泛泛地说就是相互干扰；相干解调是指利用乘法器，输入一路与载频相干（同频同相）的参考信号与载频相乘。

这里的解调是运用的非相干解调：

我们将上述的调制信号分别乘以coswt,-sinwt,就可以得到相应的IQ路的数据，然后分别通过一个低通滤波器，滤除高频分量，就可以的得到以基带信号为频率的I路Q路载波信号，然后两个再相比，就可以得到其正切值，再求反正切就可以了的到其角度，角度就是我们的基带信号。

这里的微分在离散域其实就是我们的相连的数据求差分。所以我们直接就可以求得基带信号m.

FM调制matlab实现。

说明：因为后续会用硬件FPGA 实现，所以matlab 实现实现算法也会像尽可能的逼近其硬件。

思路：

产生基本的信号：基带信号m(n),载波cos(wt),sin(wt)的实现

我们这里实现信号发生器是用的查表法，即FPGA，MATLAB 产生DDS,这里就不说明其方法了，相关信息可以看我发的其他博客，里面有详细的教程，以及原理介绍。

下面就是FM调制解调代码：

复制代码

clc;
clear all;
%*****************************************
fs = 16e6; %采样率 载波的采样率
fc = 1e6; %载波中心频率
df = 75e3; %最大频偏
fm = 16e3; %音频的采样率
kf = (df*2^32/fs)/32767;
AC= 1024;%幅度
%由于我们的载波的采样率是我们的音频的采样率的1000倍，但是我们在进行FM调制时
%会用到两个相乘，所以相当于一个语音信号的点，要被采样1000次，并且要求数据长度一致
%所以我们可以将语音信号的每一个样点复制1000次
%read pcm file 
%*********************
%read pcm file
file= fopen('v3edu.pcm','r');
m=fread(file,'int16');

m_len=length(m)*fs/fm;%按照载波采样率的长度计算
figure(1);
plot(m);title('原始的语音信号');
%*****************************
%generate 1Mhz
w=fc*2^32/fs; %频率控制字，32位的相位累加器， f_out  = fs*w/2^N :
% 载波的ROM
n=0:1/1024:1023/1024;
s_rom=sin(2*pi*n);
c_rom=cos(2*pi*n);
%plot(s_rom);
%接下来产生载波地址，
w_r=0;%相位累加器
rrom_addr=0;%
cw_sin=zeros(1,m_len);%这个相当于是我们的查找表，一个周期的正弦波
cw_cos=zeros(1,m_len);
for i=1:m_len
    w_r = w_r + w;
    if(w_r > 2^32) % 做32位累加器的溢出判断
        w_r= w_r - 2^32;
    end
    rrom_addr=round(w_r/2^22);%读查找表的地址
    if rrom_addr == 0%matlab计数时从1开始计数
        rrom_addr =1;
    end
    cw_sin(i)=s_rom(rrom_addr);%%%%%%载波完成
    cw_cos(i)=c_rom(rrom_addr);%载波
end
figure(2);
subplot(2,1,1);plot(cw_sin(1:3000),'r');title('sin载波信号');

subplot(2,1,2);plot(cw_cos(1:3000),'b');title('cos载波信号');
%给音频数据升采样，与载波长度一致，后面要相乘，所以要保证信号的长度是一样的
m_t=zeros(1,m_len);
for i=1:length(m)
    for j=1:fix(fs/fm)
        m_t((i-1)*fix(fs/fm)+j)=m(i);%同一个16k的采样点复制1000次这样就是16M采样点
    end
end

%*************************
%调制 sin 和cos 的信号
w_r=0;%相位累加器
rrom_addr=0;%
dac_i=zeros(1,m_len);%cos，I路信号 
dac_q=zeros(1,m_len);%sin  Q路信号
for i=1:m_len
    w_r = w_r + kf*m_t(i);
    if(w_r > 2^32) % 做32位累加器的溢出判断
        w_r= w_r - 2^32;
    elseif(w_r <0)
        w_r = w_r + 2^32; % 负的溢出时
    end
    
    rrom_addr=round(w_r/2^22);%读查找表的地址
    if rrom_addr == 0
        rrom_addr =1;
    end
    dac_q(i)=AC*s_rom(rrom_addr);
    dac_i(i)=AC*c_rom(rrom_addr);%载波
end

%****************************
%正交调制
s_t=zeros(1,m_len);
for i=1:m_len
    s_t(i) = dac_i(i)*cw_cos(i) + dac_q(i)*cw_sin(i)*(-1);
end
figure(3);
plot(s_t);title('调制完成信号');
%%%%%%%%%%%%调制结束%%%%%%%%%%%%%%%%%%%%%%%%%%%

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
clc;
clear all;
%*****************************************
fs = 16e6; %采样率 载波的采样率
fc = 1e6; %载波中心频率
df = 75e3; %最大频偏
fm = 16e3; %音频的采样率
kf = (df*2^32/fs)/32767;
AC= 1024;%幅度
%由于我们的载波的采样率是我们的音频的采样率的1000倍，但是我们在进行FM调制时
%会用到两个相乘，所以相当于一个语音信号的点，要被采样1000次，并且要求数据长度一致
%所以我们可以将语音信号的每一个样点复制1000次
%read pcm file 
%*********************
%read pcm file
file= fopen('v3edu.pcm','r');
m=fread(file,'int16');

m_len=length(m)*fs/fm;%按照载波采样率的长度计算
figure(1);
plot(m);title('原始的语音信号');
%*****************************
%generate 1Mhz
w=fc*2^32/fs; %频率控制字，32位的相位累加器， f_out  = fs*w/2^N :
% 载波的ROM
n=0:1/1024:1023/1024;
s_rom=sin(2*pi*n);
c_rom=cos(2*pi*n);
%plot(s_rom);
%接下来产生载波地址，
w_r=0;%相位累加器
rrom_addr=0;%
cw_sin=zeros(1,m_len);%这个相当于是我们的查找表，一个周期的正弦波
cw_cos=zeros(1,m_len);
for i=1:m_len
    w_r = w_r + w;
    if(w_r > 2^32) % 做32位累加器的溢出判断
        w_r= w_r - 2^32;
    end
    rrom_addr=round(w_r/2^22);%读查找表的地址
    if rrom_addr == 0%matlab计数时从1开始计数
        rrom_addr =1;
    end
    cw_sin(i)=s_rom(rrom_addr);%%%%%%载波完成
    cw_cos(i)=c_rom(rrom_addr);%载波
end
figure(2);
subplot(2,1,1);plot(cw_sin(1:3000),'r');title('sin载波信号');

subplot(2,1,2);plot(cw_cos(1:3000),'b');title('cos载波信号');
%给音频数据升采样，与载波长度一致，后面要相乘，所以要保证信号的长度是一样的
m_t=zeros(1,m_len);
for i=1:length(m)
    for j=1:fix(fs/fm)
        m_t((i-1)*fix(fs/fm)+j)=m(i);%同一个16k的采样点复制1000次这样就是16M采样点
    end
end

%*************************
%调制 sin 和cos 的信号
w_r=0;%相位累加器
rrom_addr=0;%
dac_i=zeros(1,m_len);%cos，I路信号 
dac_q=zeros(1,m_len);%sin  Q路信号
for i=1:m_len
    w_r = w_r + kf*m_t(i);
    if(w_r > 2^32) % 做32位累加器的溢出判断
        w_r= w_r - 2^32;
    elseif(w_r <0)
        w_r = w_r + 2^32; % 负的溢出时
    end
    
    rrom_addr=round(w_r/2^22);%读查找表的地址
    if rrom_addr == 0
        rrom_addr =1;
    end
    dac_q(i)=AC*s_rom(rrom_addr);
    dac_i(i)=AC*c_rom(rrom_addr);%载波
end

%****************************
%正交调制
s_t=zeros(1,m_len);
for i=1:m_len
    s_t(i) = dac_i(i)*cw_cos(i) + dac_q(i)*cw_sin(i)*(-1);
end
figure(3);
plot(s_t);title('调制完成信号');
%%%%%%%%%%%%调制结束%%%%%%%%%%%%%%%%%%%%%%%%%%%

MATLAB 滤波器的实现：

我们实现一个简单的低通滤波器，只需要将我们我们需要的滤波器参数输入Filter Solutions中，然后复制器抽头系数，然后与我们的输入数据相卷积就可以了，下面的代码有注释.

复制代码

%%%%%%%%%%%%正交解调%%%%%%%%%%%%%%%%
i_data = zeros(1,m_len);
q_data = zeros(1,m_len);
for i = 1 : m_len
	i_data(i) = s_t(i) * cw_cos(i);
	q_data(i) = (-1)* s_t(i) * cw_sin(i);
end
%%%%%%%%%%%%%%进行低通滤波%%%%%%%%%%%%%%%%
%滤波器抽抽头系数，
NUM = [-7.894e-05, -2.483e-04, -4.516e-04, -7.089e-04, -1.035e-03, -1.435e-03, -1.903e-03, -2.417e-03, -2.938e-03, -3.415e-03, -3.778e-03, -3.948e-03, -3.837e-03, -3.357e-03, -2.42e-03, -9.524e-04, 1.106e-03, 3.791e-03, 7.113e-03, 1.105e-02, 1.554e-02, 2.051e-02, 2.583e-02, 3.135e-02, 3.691e-02, 4.233e-02, 4.742e-02, 5.199e-02, 5.588e-02, 5.894e-02, 6.105e-02, 6.213e-02, 6.213e-02, 6.105e-02, 5.894e-02, 5.588e-02, 5.199e-02, 4.742e-02, 4.233e-02, 3.691e-02, 3.135e-02, 2.583e-02, 2.051e-02, 1.554e-02, 1.105e-02, 7.113e-03, 3.791e-03, 1.106e-03, -9.524e-04, -2.42e-03, -3.357e-03, -3.837e-03, -3.948e-03, -3.778e-03, -3.415e-03, -2.938e-03, -2.417e-03, -1.903e-03, -1.435e-03, -1.035e-03, -7.089e-04, -4.516e-04, -2.483e-04, -7.894e-05];
adc_i = conv(i_data,NUM); %低通滤波
adc_q = conv(q_data,NUM);%%%%
%%%%%%%%%%%%%FM解调%%%%%%%%%%%%%%%%%%%%%%
%重新检测数据长度
c_len = length(adc_i);
cr = zeros(1,c_len);
cj = zeros(1,c_len);
for i = 2: c_len %必须从2开始，因为有I-1 
	cr(i) = adc_i(i)*adc_i(i) +adc_q(i)*adc_q(i);
	cj(i)  = adc_i(i-1) * adc_q(i) -  adc_i(i)* adc_q(i-1);%I(n-1)*Q(n) -I(n)*Q(n-1)
end
angle = zeros(1,c_len);
for i = 1:c_len
	if cr(i) == 0
		angle(i) = 0;
	else 
		angle(i) = (cj(i)/cr(i))*2 ^16; %FPGA 中没有小数，要进行量化
	end
end
%%%%%%%%%%%%%降采样%%%%%%%%%%%%%%%%%%%%%%%%%
%降采样，以均值方式，16M，降到16k
d_len = length(angle);
sum = 0;
cnt = 0;
demout = zeros(1,fix(d_len*(fm/fs)));
d_cnt = 1;
for i= 1:d_len
	sum = sum + angle(i);
	if cnt == fs/fm %完成1000次的累加
		demout(d_cnt) =  sum /cnt;
		d_cnt = d_cnt + 1;
		cnt = 0;
		sum= 0;
	end
	cnt = cnt + 1;
end
figure(4);
plot(demout);title('解调后的信号');

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
%%%%%%%%%%%%正交解调%%%%%%%%%%%%%%%%
i_data = zeros(1,m_len);
q_data = zeros(1,m_len);
for i = 1 : m_len
	i_data(i) = s_t(i) * cw_cos(i);
	q_data(i) = (-1)* s_t(i) * cw_sin(i);
end
%%%%%%%%%%%%%%进行低通滤波%%%%%%%%%%%%%%%%
%滤波器抽抽头系数，
NUM = [-7.894e-05, -2.483e-04, -4.516e-04, -7.089e-04, -1.035e-03, -1.435e-03, -1.903e-03, -2.417e-03, -2.938e-03, -3.415e-03, -3.778e-03, -3.948e-03, -3.837e-03, -3.357e-03, -2.42e-03, -9.524e-04, 1.106e-03, 3.791e-03, 7.113e-03, 1.105e-02, 1.554e-02, 2.051e-02, 2.583e-02, 3.135e-02, 3.691e-02, 4.233e-02, 4.742e-02, 5.199e-02, 5.588e-02, 5.894e-02, 6.105e-02, 6.213e-02, 6.213e-02, 6.105e-02, 5.894e-02, 5.588e-02, 5.199e-02, 4.742e-02, 4.233e-02, 3.691e-02, 3.135e-02, 2.583e-02, 2.051e-02, 1.554e-02, 1.105e-02, 7.113e-03, 3.791e-03, 1.106e-03, -9.524e-04, -2.42e-03, -3.357e-03, -3.837e-03, -3.948e-03, -3.778e-03, -3.415e-03, -2.938e-03, -2.417e-03, -1.903e-03, -1.435e-03, -1.035e-03, -7.089e-04, -4.516e-04, -2.483e-04, -7.894e-05];
adc_i = conv(i_data,NUM); %低通滤波
adc_q = conv(q_data,NUM);%%%%
%%%%%%%%%%%%%FM解调%%%%%%%%%%%%%%%%%%%%%%
%重新检测数据长度
c_len = length(adc_i);
cr = zeros(1,c_len);
cj = zeros(1,c_len);
for i = 2: c_len %必须从2开始，因为有I-1 
	cr(i) = adc_i(i)*adc_i(i) +adc_q(i)*adc_q(i);
	cj(i)  = adc_i(i-1) * adc_q(i) -  adc_i(i)* adc_q(i-1);%I(n-1)*Q(n) -I(n)*Q(n-1)
end
angle = zeros(1,c_len);
for i = 1:c_len
	if cr(i) == 0
		angle(i) = 0;
	else 
		angle(i) = (cj(i)/cr(i))*2 ^16; %FPGA 中没有小数，要进行量化
	end
end
%%%%%%%%%%%%%降采样%%%%%%%%%%%%%%%%%%%%%%%%%
%降采样，以均值方式，16M，降到16k
d_len = length(angle);
sum = 0;
cnt = 0;
demout = zeros(1,fix(d_len*(fm/fs)));
d_cnt = 1;
for i= 1:d_len
	sum = sum + angle(i);
	if cnt == fs/fm %完成1000次的累加
		demout(d_cnt) =  sum /cnt;
		d_cnt = d_cnt + 1;
		cnt = 0;
		sum= 0;
	end
	cnt = cnt + 1;
end
figure(4);
plot(demout);title('解调后的信号');