function [msk_sig sig_n] = mskmod(sig,fs,fb,fc,snr)
%**********初始化**********
if nargin < 5
snr = 10; %信噪比，以dB为单位
end
if nargin < 4
fc = 2; %载频
end;
if nargin < 3
fb = 2; %码元速率
end;
if nargin < 2
fs = 200; %采样频率
end;
if nargin < 1
sig = randint(1,10,2); %待调制信号
end;
ls = length(sig);
a = mod(ls,2);
if a ~= 0
error('Please insert a vector divisible for 2');
end
m = fs/fb; %每个码元上的采样频率
T = length(sig)/fb; %仿真总时间，信号点个数*每个数保持的时间（即码元长度）
dt = 1/fs; %两个采样点之间的时间间隔
t = 0:dt:T-dt;
sig1 = repmat(sig,m,1);
sig2 = reshape(sig1,1,[]); %此两步用于生成二进制方波，每个二进制数保持m次
I=zeros(1,ls);
Q=zeros(1,ls);
phi=zeros(1,ls); %初始化相位常数φ
for k=2:ls
phi(k)=phi(k-1)+2*(sig(k-1)-sig(k))*((k-1)*pi/2);
end
%*********定义方法************
I=cos(phi);
Q=(sig*2-1).*I; %Q=ak*I,其中sig*2-1实现单极性转变到双极性
I1 = repmat(I,m,1);
I2 = reshape(I1,1,[]);
Q1 = repmat(Q,m,1);
Q2 = reshape(Q1,1,[]);
msk_sig = I2.*cos(pi*t*fb/2).*cos(2*pi*fc*t)-Q2.*sin(pi*t*fb/2).*sin(2*pi*fc*t);
%**********加高斯白噪声**********
sig_n = awgn(msk_sig,snr,'measured');
%**********画图*****************
subplot(2,1,1);
plot(sig2,'LineWidth',1.5);
grid on;
title('二进制信号');
axis([0 m*length(sig) -2.5 2.5]);
subplot(2,1,2);
plot(t,sig2,t,msk_sig,'LineWidth',1.5);
grid on;
title('MSK调制信号');
axis([0 T -2.5 2.5]);
评论0