mskgmsk.m

利用matlab仿真GMSK的源程序 按照移动通信的原理编写的程序

%绘制调制波形10110001
clear all;
Ts=1/16000;                         %基带信号周期为1/16000s,即为16KHz
Tb=1/32000;                         %输入信号周期为Ts/2=1/32000s,即32KHz
BbTb=0.5;                           %取BbTb为0.3
Bb=BbTb/Tb;                         %3dB带宽
Fc=32000;                           %载波频率为32KHz
F_sample=64;                        %每载波采样64个点
B_num=8;                            %基带信号为8个码元
B_sample=F_sample*Fc*Tb             %每基带码元采样点数B_sample=Tb/Dt
Dt=1/Fc/F_sample;                   %采样间隔
t=0:Dt:B_num*Tb-Dt;                 %仿真时间
T=Dt*length(t);                     %仿真时间值
%Data=sign(randn(1,B_num));         %随机产生100个基带信号
Ak=[1 0 1 1 0 0 0 1];               %产生8个基带信号
Ak=2*Ak-1;
gt=ones(1,B_sample);         %每码元对应的载波信号
Akk=sigexpand(Ak,B_sample);   %码元扩展
temp=conv(Akk,gt);                 %码元扩展
Akk=temp(1:length(Akk));            %码元扩展

Tbt=0:Dt:Tb-Dt;                 %每Tb时间
ThetaL=0;
Theta=zeros(1,length(Akk));
for k=1:B_num
    if k==1
        ThetaL=0;
    else
        ThetaL=Theta((k-1)*B_sample);
    end
    for j=1:B_sample
        Theta((k-1)*B_sample+j)=pi*Ak(k)*Tbt(j)/2/Tb+ThetaL;
    end    
end;

S_Msk=cos(2*pi*Fc*t+Theta);
[f F_Msk]=T2F(t,S_Msk);

tt=-2.5*Tb:Dt:2.5*Tb-Dt;   
%g(t)=Q[2*pi*Bb*(t-Tb/2)/sqrt(log(2))]-Q[2*pi*Bb*(t+Tb/2)/sqrt(log(2))];
%Q(t)=erfc(t/sqrt(2))/2;
gausst=erfc(2*pi*Bb*(tt-Tb/2)/sqrt(log(2))/sqrt(2))/2-erfc(2*pi*Bb*(tt+Tb/2)/sqrt(log(2))/sqrt(2))/2;    

J_g=zeros(1,length(gausst); %使J_g 的长度和Gausst的一样
for i=1:length(gausst)
    if i==1 
        J_g(i)=gausst(i)*Dt;
    else
        J_g(i)=J_g(i-1)+gausst(i)*Dt;
    end;
end;
J_g=J_g/2/Tb;

Alpha=zeros(1,length(Akk));
k=1;
L=0;
for j=1:B_sample
    J_Alpha=Ak(k+2)*J_g(j);
    Alpha((k-1)*B_sample+j)=pi*J_Alpha+L*pi/2;
end; 

k=2;
L=0;
for j=1:B_sample
    J_Alpha=Ak(k+2)*J_g(j)+Ak(k+1)*J_g(j+B_sample);
    Alpha((k-1)*B_sample+j)=pi*J_Alpha+L*pi/2;
end;  

k=3;
L=0;
for j=1:B_sample
    J_Alpha=Ak(k+2)*J_g(j)+Ak(k+1)*J_g(j+B_sample)+Ak(k)*J_g(j+2*B_sample);
    Alpha((k-1)*B_sample+j)=pi*J_Alpha+L*pi/2;
end;  

k=4;
L=0;
for j=1:B_sample
    J_Alpha=Ak(k+2)*J_g(j)+Ak(k+1)*J_g(j+B_sample)+Ak(k)*J_g(j+2*B_sample)+Ak(k-1)*J_g(j+3*B_sample);
    Alpha((k-1)*B_sample+j)=pi*J_Alpha+L*pi/2;
end;

L=0;
for k=5:B_num-2
    if k==5
        L=0;
    else
        L=L+Ak(k-3);
    end;
    for j=1:B_sample
        J_Alpha=Ak(k+2)*J_g(j)+Ak(k+1)*J_g(j+B_sample)+Ak(k)*J_g(j+2*B_sample)+Ak(k-1)*J_g(j+3*B_sample)+Ak(k-2)*J_g(j+4*B_sample);
        Alpha((k-1)*B_sample+j)=pi*J_Alpha+mod(L,4)*pi/2;
    end;    
end;

%B_num-1;
k=B_num-1;
L=L+Ak(k-3);
for j=1:B_sample
    J_Alpha=Ak(k+1)*J_g(j+B_sample)+Ak(k)*J_g(j+2*B_sample)+Ak(k-1)*J_g(j+3*B_sample)+Ak(k-2)*J_g(j+4*B_sample);
    Alpha((k-1)*B_sample+j)=pi*J_Alpha+mod(L,4)*pi/2;
end;  
%B_num;
k=B_num;
L=L+Ak(k-3);
for j=1:B_sample
    J_Alpha=Ak(k)*J_g(j+2*B_sample)+Ak(k-1)*J_g(j+3*B_sample)+Ak(k-2)*J_g(j+4*B_sample);
    Alpha((k-1)*B_sample+j)=pi*J_Alpha+mod(L,4)*pi/2;
end;  


S_Gmsk=cos(2*pi*Fc*t+Alpha);
[f F_Gmsk]=T2F(t,S_Gmsk);

figure(1);
subplot(2,1,1);
plot(tt/Tb,gausst);
axis([-2.5 2.5 0 1]);
subplot(2,1,2);
plot(tt/Tb,J_g);
axis([-2.5 2.5 0 0.5]);


Show_Num=8;                     %显示码元数
Show_Time=Show_Num*Tb;          %显示码元数

figure(2);
subplot(421)
plot(t/Tb,Akk);
axis([0 Show_Num -1.5 1.5]);
title('基带波形');

subplot(423)
plot(t/Tb,Theta*2/pi);
axis([0 Show_Num min(Theta*2/pi)-1 max(Theta*2/pi)+1]);
title('相位');

subplot(425)
plot(t/Tb,S_Msk);
axis([0 Show_Num -1.5 1.5]);
title('MSK波形');

subplot(427)
plot(f/Fc,10*log10(abs(F_Msk).^2/length(F_Msk)));
axis([-15 15 -200 -100]);
xlabel('f');
ylabel('MSK功率谱密度(dB/Hz)');
title('MSK功率谱密度');

subplot(422)
plot(t/Tb,Akk);
axis([0 Show_Num -1.5 1.5]);
title('基带波形');

subplot(424)
plot(t/Tb,Alpha*2/pi);
axis([0 Show_Num min(Alpha*2/pi)-1 max(Alpha*2/pi)+1]);
title('相位波形');

subplot(426)
plot(t/Tb,S_Gmsk);
axis([0 Show_Num -1.5 1.5]);
title('GMSK波形');

subplot(428)
plot(f/Fc,10*log10(abs(F_Gmsk).^2/length(F_Gmsk)));
axis([-15 15 -200 -100]);
xlabel('f');
title('GMSK功率谱密度');