MATLAB通信系统仿真（一）模拟调制与解调

模拟调制与解调

---

---

前言

主要内容来自参考资料[2]，在学习记录的同时勘正了书中代码的部分错误。

---

一、拟合

常用函数：polyfit()，polyval()，功能如下图所示：

二、调制与解调仿真

代码如下（示例）：

SSB（单边带调制）

希尔伯特变换

$$\hat{m}(t)=m(t)/\pi t\hat{M}(f)=-j\times sgn(f)M(f)$$其中$\hat{m}(t)$为$m(t)$的希尔伯特变换，而$\hat{M}(f)$，$M(f)$为两者的傅里叶变换

假设$m(t)$为消息信号，则对其进行单边带调制，得到调制信号$u(t)$为
$$u(t)=A_{c}m(t)cos(2\pi f_{c}t)/2\mp A_c\hat{m}(t)cos(2\pi f_{c}t)/2$$使用MATLAB进行仿真。对正弦信号 y = sin(t) 进行希尔伯特变换，并绘图。

clear all;close all
t = 0:0.1:30;
y = sin(t);
s_y = hilbert(y);
plot(t,y,'linewidth',2);hold on
plot(t,real(s_y),'ro',t,imag(s_y),'bs');
legend('原信号','希尔伯特变换器实部（仍然是原信号）','希尔伯特变换器虚部')

绘图结果如下

SSB调制

基带信号是在150Hz~400Hz内，幅度随频率逐渐递减的音频信号，对其进行SSB调制,仿真时间为1s。载波频率为1000Hz，采样频率为10000Hz。

示例MATLAB代码：

% 基带信号是在150Hz~400Hz内，幅度随频率逐渐递减的音频信号，对其进行SSB调制,仿真时间为1s
clear all;close all;
fs = 10000;
fc = 1000;
t = 1/fs:1/fs:1;%仿真时间
f_v = 0:length(t)-1;
m(fs)=0;%创建长度为fs的消息信号
for f = 150:400
    m = m + 0.01*sin(2*pi*f*t)*(400-f);%m(t)包含许多频率分量，且频率分量的幅度随频率递减
end
m_shift = imag(hilbert(m));%基带信号的希尔伯特变换
carriercos = cos(2*pi*fc*t);%载波
carriersin = sin(2*pi*fc*t);%载波的正交变换
ssb_u = m.*carriercos - m_shift.*carriersin; %上边带SSB
ssb_l = m.*carriercos + m_shift.*carriersin; %上边带SSB

figure;
subplot(421);plot(t,carriercos,t,carriersin,':m');axis([0 0.01 -1 1]);title('载波及其正交变换')%载波
subplot(422);plot(f_v,abs(fft(carriercos)));axis([0 2000 -500 12000]);title('载波频谱图')      %载波频谱
subplot(423);plot(t,m);axis([0 0.01 -500 500]);title('基带信号波形图')                       	%基带信号
subplot(424);plot(f_v,abs(fft(m)));axis([0 2000 -500 12000]);title('基带信号频谱图')        	%基带信号频谱图
subplot(425);plot(t,ssb_u);axis([0 0.01 -500 500]);title('%SSB上边带波形图')                   %SSB上边带波形图
subplot(426);plot(f_v,abs(fft(ssb_u)));axis([0 2000 -500 12000]);title('SSB上边带频谱图')      %SSB上边带频谱图
subplot(427);plot(t,ssb_l);axis([0 0.01 -500 500]);title('SSB下边带波形图')                    %SSB下边带波形图
subplot(428);plot(f_v,abs(fft(ssb_l)));axis([0 2000 -500 12000]);title('SSB下边带频谱图')      %SSB下边带频谱图
suptitle('波形图与频谱图')

仿真结果：

SSB解调

在实际应用中，一般采用高稳定度的晶体振荡器或频率合成器来产生本地解调载波，无需像DSB解调那样用PLL恢复载波，大大降低了单边带接收机的技术复杂度与成本。

在本例中，将采用相干解调法，即产生本地振荡波形与接收信号相乘，再经过低通滤波，恢复原始信号，但由于本地振荡与载波不严格同频同相，此法存在一定的频率偏移和相位偏移。

示例代码：

out = ssb_u.*carriercos;%相干解调
[a,b] = butter(4,500/(fs/2));    %设计4阶低通滤波器，且截止频率为500Hz
out1 = filter(a,b,out);       %滤波输出

figure;
subplot(323);plot(t,out);axis([0 0.01 -500 500]);title('解调信号波形图')                       	
subplot(324);plot(f_v,abs(fft(out)));axis([0 2000 -500 12000]);title('解调信号频谱图')        
subplot(325);plot(t,out1);axis([0 0.01 -500 500]);title('滤波信号波形图')                       
subplot(326);plot(f_v,abs(fft(out1)));axis([0 2000 -500 12000]);title('滤波信号频谱图')        
subplot(321);plot(t,ssb_u);axis([0 0.01 -500 500]);title('%SSB上边带波形图')                   %SSB上边带波形图
subplot(322);plot(f_v,abs(fft(ssb_u)));axis([0 2000 -500 12000]);title('SSB上边带频谱图')      %SSB上边带频谱图
suptitle('波形图与频谱图')

仿真结果图：

角度调制

频率调制

示例代码

%% 频率调制 FM
clear all;close all;
t0 = 0.15;%信号保持时间
ts = 0.0005;%采样时间
fc = 200;%载波频率
kf = 50;%调制系数
fs = 1/ts;%采样频率
df = 0.25;%频率分辨率
t_v = [0:ts:t0];%时间序列
m = [ones(1,t0/(3*ts)),-2*ones(1,t0/(3*ts)),zeros(1,t0/(3*ts)+1)];%消息信号m(t)，t/ts是序列的数目。注意端点，要+1
m_int(length(m)) = 0;
%以下两种方法都可实现对m(t)的积分。法二更直观，一定要乘以ts！
for i = 1:length(m)-1
    m_int(i+1) = m_int(i) + m(i)*ts;
end
% for i = 1:length(m)
%     m_int(i) = ts*sum(m(1:i));
% end
[M,m,df1] = fftseq(m,ts,df);%傅里叶变换
M = M/fs;
f = [0:df1:df1*(length(m)-1)]-fs/2;%画[-fs/2,fs/2]的图像
u = cos(2*pi*fc*t_v + 2*pi*kf*m_int);%调制信号
figure
plot(t_v,u)
[U,u,df1] = fftseq(u,ts,df);%傅里叶变换
U = U/fs;
%画图
figure;
subplot(221);plot(t_v,m(1:length(t_v)));%经过fftseq函数后，m的长度与t_v的长度不一定相等，但有效长度相等
axis([0 0.15 -2.1 2.1]);
xlabel('时间');title('消息信号')
subplot(222);plot(t_v,u(1:length(t_v)))
axis([0 0.15 -2.1 2.1]);
xlabel('时间');title('调频信号')

subplot(223);plot(f,abs(fftshift(M)));
xlabel('频率');title('消息信号频谱')
subplot(224);plot(f,abs(fftshift(U)));
xlabel('频率');title('调制信号频谱')

仿真结果图：

相位调制

示例代码

clear all;close all;
t0 = 0.25;%信号保持时间
ts = 0.0005;%采样时间
fc = 200;%载波频率
fs = 1/ts;%采样频率
df = 0.25;%频率分辨率
t_v = [0:ts:t0];%时间序列
m(length(t_v)) = 0;
lent = length(t_v)-1;
for i = 1:(t0/ts/4)
    m(i) = i;
end
for i = (lent/4+1):3*(lent/4)
    m(i) = m(lent/4) - i+lent/4;
end
for i = (3*lent/4+1):lent+1
    m(i) = m(3*lent/4) + i-lent*3/4
end
m = m/50;

[M,m,df1] = fftseq(m,ts,df);%傅里叶变换
M = M/fs;
f = [0:df1:df1*(length(m)-1)]-fs/2;%画[-fs/2,fs/2]的图像

u = cos(2*pi*fc*t_v + m(1:length(t_v)));%调制信号
[U,u,df1] = fftseq(u,ts,df);%傅里叶变换
U = U/fs;
%画图
figure;
subplot(221);plot(t_v,m(1:length(t_v)));%经过fftseq函数后，m的长度与t_v的长度不一定相等，但有效长度相等
axis([0 0.25 -3 3]);
xlabel('时间');title('消息信号')
subplot(222);plot(t_v,u(1:length(t_v)))
axis([0 0.15 -2.1 2.1]);
xlabel('时间');title('调频信号')

subplot(223);plot(f,abs(fftshift(M)));
xlabel('频率');title('消息信号频谱')
subplot(224);plot(f,abs(fftshift(U)));
xlabel('频率');title('调制信号频谱')

仿真结果图：

相关函数

fft()和fftshift()的区别

M = fft(m)是对消息信号 m(t) 进行快速FFT变换，得到取值范围为 (0,fs) 的函数 M(f) ，如果想画出取值范围为 (-fs/2,fs/2) 的函数时，需要用到 fftshift() 函数。

需要注意的是，使用时直接 fftshift(M) 即可，该函数只是实现对函数 M(f) 的移位，并未对其进行FFT变换。

详情见博客： https://www.cnblogs.com/limanjihe/p/10014142.html

参考资料

[1] ttps://blog.csdn.net/u014630987/article/details/70156489
[2] 《MATLAB通信系统建模与仿真（第2版）》 邓奋发 编著