MSK、GMSK调制及相干解调Matlab仿真

MSK、GMSK调制及相干解调Matlab仿真实验目的：

1.掌握MSK调制、相干解调原理及特性；

2.了解MSK调制与GMSK调制的差别。

实验内容：

1.编写MATLAB程序仿真MSK调制及相干解调；

2.观察I、Q两路基带信号的特征及与输入NRZ码的关系；

3.观察I、Q调制解调过程中信号的变化；

4.对程序做修改，进行GMSK调制及解调仿真；

5.分析仿真中观察的数据，撰写实验报告。

仿真代码：

clear all

clc

global dt df t f N

close all

pi=3.1415926;

fc=5;

N=2^8;

L=8;

M=N/L;

Rb=2;

Tb=1/Rb;

dt=Tb/L;

df=1/(dt*N);

T=N*dt;

B=N*df/2;

t=[-T/2+dt/2:dt:T/2];

f=[-B+df/2:df:B];

EP=zeros(size(f));

EPg=zeros(size(f));

for ii=1:10;

for j=1:50;

b=sign(randn(1,M));

for i=1:L,s(i+[0:M-1]*L)=b;end

P=t2f(s);

P=P.*conj(P)/T;

EP=(EP*(j-1)+P)/j;

end

Ps=10*log10(EP+eps);

Bb=Tb/0.3;

alpha=sqrt(logm(2)/2/Bb^2);

H=exp(-alpha^2*f.^2);

a(1)=b(1);

for i=M:-1:2,a(i)=b(i)*b(i-1);end

for i=1:L,sa(i+[0:M-1]*L)=a;end

send=real(f2t(t2f(s).*H));

It=zeros(size(t));

for k=0:2*L:N-1;

kk=1:2:2*L;

kkk=1:L;

It(k+kk)=send(k+kkk+L);

It(k+kk+1)=send(k+kkk+L);

end

for k=N:-1:L+1,It(k)=It(k-L);end

Qt=zeros(size(t));

for k=0:2*L:N-1;

kk=1:2:2*L;

kkk=1:L;

Qt(k+kk)=send(k+kkk);

Qt(k+kk+1)=send(k+kkk);

end

subplot(2,1,1);

stem(b);

title('原始');

%x = input('xxx');

%Itt

Itt=It.*cos(pi*t/2/Tb);

%Qtt

Qtt=Qt.*sin(pi*t/2/Tb);

%GMSK 时域波形

gmsk=Itt.*cos(2*pi*fc*t)-Qtt.*sin(2*pi*fc*t); %GMSK 功率谱

PP=t2f(gmsk);

Pa=PP.*conj(PP)/T;

EPg=(EPg*(ii-1)+Pa)/ii;

end

Pgmsk=10*log10(EPg+eps);

%接收端

r=gmsk;

%接收端的低通滤波器，带宽为Rb

LPF=zeros(size(f));

ai=(B-Rb)/2/B*size(f);

aj=(B+Rb)/2/B*size(f);

for k=(ai(1,2):aj(1,2)),LPF(k)=1;end

%接收端上支路LPF的输出，与Itt相似（图九）RI=r.*cos(2*pi*fc*t);

RI=real(f2t(t2f(RI).*LPF));

RQ=-r.*sin(2*pi*fc*t);

RQ=real(f2t(t2f(RQ).*LPF));

%取样

RIt=RI(2*L:2*L:N);

RQt=RQ(L:2*L:N);

%码型串并转换

Rt=zeros(1,M);

Rt(2:2:M)=RIt(1:M/2);

Rt(1:2:M-1)=RQt(1:M/2);

%判决

Rt=sign(Rt);

clear j;

d(1)=j;

for i=2:M,d(i)=d(i-1).*j;end

e=Rt.*d;

for i=1:2:M,e(i)=imag(e(i));end

f=b-e;

for i=1:L,sy(i+[0:M-1]*L)=e;end

subplot(2,1,2);

stem(sy);

f2t

function x=f2t(X)

global dt df t f T N

X=[X(N/2+1:N),X(1:N/2)];

x=ifft(X)/dt;

t2f

function X=t2f(x)

global dt N t f T

H=fft(x);

X=[H(N/2+1:N),H(1:N/2)]*dt;