MSK调制解调的Matlab仿真

整理了一下MSK调制解调的Matlab仿真仅供大家参考，

MSK调制解调的基本原理请参考本人的记录：

MSK调制解调原理_msk调制及解调实验电路原理-CSDN博客

上面文章中很多需要自己理解的东西无法写的很详细，需要自己通过仿真和实验来体会。这里将自己整理的matlab仿真提供给大家学习参考。大家可以结合仿真代码和MSK调制解调的原理框图来学习。

代码纯手敲，注释详细，如果对大家有帮助，麻烦大家点赞、评论和收藏。O(∩_∩)O

% msk调制解调matlab仿真:验证差分编码，和MSK调制解调
clear;close all;clc;

N   =  32;                         % 一个码元中bit数:32bit
Tb  =  0.2e-6;                     % 比特周期，一个码元周期为0.2us
fc  =  20e6;                       % 载波频率
fs  =  100e6;                      % 采样频率
snr =  8;                          % 信躁比

num = fs*Tb*N;                     % 采样总点数
t=linspace(0,N*Tb,num);
% ********************** 生成测试比特流 ********************************%
%bitstream=[1,0,1,0,0,1,1,0,1];
%bitstream=[1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0];
bitstream=randi([0,1],1,N);       % 随机生成比特流


% ************************************************************************%
%%  1、信号调制 
% ************************************************************************%

% ************************ 1.1双极性编码 **********************************%

ak=2*bitstream-1;                  % 产生双极性比特流 -1,1

% ************************* 1.2差分编码 *********************************%
bk(1)=ak(1);
for i=2:N
    bk(i)=ak(i)*bk(i-1); %#ok<*SAGROW> %对应bk
end

% ************************* 1.3串并转换 *********************************%
%奇数进I路、偶数进Q路
I=[];Q=[];
for i=1:N
    if mod(i,2)~=0
        I=[I,bk(i)];
    else
        Q=[Q,bk(i)];
    end
end
%用绘图来比较I路和Q路比特流
%脉冲成型，认为输入的比特流是矩形的
bit_wave=[];
for i=1:N
    bit_wave=[bit_wave,bk(i)*ones(1,Tb*fs)];  %#ok<*AGROW> 
end

I_data=[];Q_data=[];
base_wave=0:1/fs:Tb-1/fs;

%************************** 1.4 脉冲成型 ***********************************%
% I路:(-T,T)之间的半圆乘I(i)是一个比特周期里面的脉冲成型
% Q路与I路一样的脉冲成型波形
% 16*(2*0.2us*2400M)=16*960=15360 

t_i = linspace (0,2*Tb,2*Tb*fs);
t_q = linspace (0,2*Tb,2*Tb*fs);
I_mod = sin(pi*t_i/(2*Tb));
Q_mod = sin(pi*t_q/(2*Tb));


for i=1:N/2
    for j = 1:2*Tb*fs
           I_data=[I_data,I(i)*I_mod(j)];       
           Q_data=[Q_data,Q(i)*Q_mod(j)];
    end
end

% % 延时一个比特周期，一个比特周期所含有的采样点个数除以2，得到一个比特周期里面点的个数
number_delay  = length(t_q)/2;
% Q路延时一个比特周期,与后面信息进行补齐;
Q_data  =  [zeros(1,number_delay),Q_data(1:length(Q_data)-number_delay)];

%得到差分编码输出的比特流,绘图与I路和Q路的信息进行比较
figure();
subplot(3,1,1);plot(t,bit_wave);xlabel('时间/us');legend('比特流');axis([0,Tb*N,-2,2]);
subplot(3,1,2);plot(I_data);xlabel('时间/us');legend('I路比特流');
subplot(3,1,3);plot(Q_data);xlabel('时间/us');legend('Q路比特流');

% ************************* 1.5 添加载波信号 ***********************************%
%载波信号
bit_t=0:1/fs:N*Tb-1/fs;%定义时间轴
%定义I路和Q路的载波信号
I_carrier = cos(2*pi*fc*bit_t);
Q_carrier = cos(2*pi*fc*bit_t + pi/2);
I_dot=I_data.*I_carrier;
Q_dot=Q_data.*Q_carrier;
MSK_signal=I_data.*I_carrier+Q_data.*Q_carrier;

figure();
t=0:1/fs:N*Tb-1/fs;
subplot(3,1,1)
plot(t,MSK_signal);legend('MSK信号：I路+Q路')
subplot(3,1,2)
plot(t,I_dot);legend('I路*Fc')
subplot(3,1,3)
plot(t,Q_dot);legend('Q路*Fc')

% ************************ 信道传输 *********************************%
% 调用awgn函数，加上噪声，得到接收的MSK信号
MSK_receive=awgn(MSK_signal,snr);
% MSK_receive = MSK_signal;
% ************************************************************************%
%%  2、信号解调 
% ************************************************************************%


%设计低通滤波器
[b,a] = butter(4,0.5,'low');

freqz(b,a,[],fs*1e-6)

%进I路时MSK接收信号与I路载波进行点乘；
I_output=MSK_receive.*I_carrier;
%调用matlab滤波函数filter把I路信号通过滤波出来得到I路滤波信号
I_filter_output=filtfilt(b,a,I_output);

%Q路与I路同理,二者一样
Q_output=MSK_receive.*Q_carrier;
Q_filter_output=filtfilt(b,a,Q_output);

figure();
subplot(3,1,1);plot(MSK_receive);xlabel('时间/us');legend('解调I路');
subplot(3,1,2);plot(I_filter_output);xlabel('时间/us');legend('解调I路');
subplot(3,1,3);plot(Q_filter_output);xlabel('时间/us');legend('解调Q路');


%抽取出I路和Q路的信息
%I路在奇数倍周期进行采样，大于0为1、小于0为-1
%Q路在偶数倍周期进行采样，大于0为1、小于0为-1
for i=1:N/2
    a = (2*i-2)*number_delay;
    b = (2*i-1)*number_delay;

    i_extract =  (a + b)/2;
    q_extract = (b + 2*i*number_delay)/2;
    x_i(i) = I_filter_output(i_extract);
    x_q(i)= Q_filter_output(q_extract);
   
    if x_i(i)>0 
        I_recover(i) = 1;
    else
        I_recover(i) = -1;
    end
    
    if  x_q(i)>0 
        Q_recover(i)=1;
    else
        Q_recover(i)=-1;
    end
end

%并/串变换得到最后输出的比特流
bit_recover=[];
for i=1:N
    if mod(i,2) ~= 0
        bit_recover=[bit_recover,I_recover((i-1)/2+1)];
    else
        bit_recover=[bit_recover,Q_recover(i/2)];
    end
end

%比特流在接收端进行差分译码
for i=1:N
    if i==1
        bit_recover_diff(i)=bit_recover(i);
    else
        bit_recover_diff(i)=bit_recover(i)*bit_recover(i-1);
    end
end

%讲双极性编码变换成单极性编码
an=[];
for i=1:N
    if bit_recover_diff(i)==1
        an(i) = 1;
    else
        an(i) = 0;
    end
end

%根据差分译码出来的数据和前面输入的比特数据进行比较，绘图看错了几个
recover_data=[];
for i=1:N
   recover_data=[recover_data,bit_recover_diff(i)*ones(1,Tb*fs)];
end
figure();
t=0:1/fs:N*Tb-1/fs;
bit_stream=[];

for i=1:N
   bit_stream=[bit_stream,ak(i)*ones(1,Tb*fs)];
end

subplot(2,1,1)
plot(t,bit_stream);legend('原始比特流'),axis([0,Tb*N,-2,2]);
subplot(2,1,2)
plot(t,recover_data);legend('恢复比特流'),axis([0,Tb*N,-2,2]);