gmskmodulationdemodulation.m

GMSK的调制解调

%图1:原始双极性不归零码
%图2:原始双极性不归零码的功率谱
%图3:高斯滤波器的传递函数
%图4:预编码后的双极性不归零码
%图5:高斯滤波后的双极性不归零码                
%图6:It的移位，将此图右移一个码元间隔(Tb=0.5)得It上支路
%图7:It(串并转换后的上支路，方波)        
%图8:Qt(串并转换后的下支路，方波) 
%图9:Itt(It与余弦相乘)
%图10:Qtt(Qt与正弦相乘)
%图11:GMSK时域信号
%图12:GMSK功率谱
%图13:接收端的低通滤波器
%图14:恢复的预编码后的双极性不归零码/接收端上支路LPF的输出，与Itt相似（图9）
%图15:接收端下支路LPF的输出,与Qtt相似（图10）
%图16:判决后的双极性不归零码
clear all
clc
global dt df t f N
close all
pi=3.1415926
fc=5    %单位:MHz
N=2^13            % 采样点数
L=64               % 每码元的采样点数
M=N/L              % 码元数
Rb=2               % 发送码元的信息速率:2Mbps
Tb=1/Rb            % 码元宽度:0.5us
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];    % 频域横坐标

figure(1);
set(1,'Position',[10,100,300,300]);   % 设定窗口位置及大小
figure(2);
set(2,'Position',[330,100,300,300]);  % 设定窗口位置及大小
figure(3);
set(3,'Position',[650,100,300,300]);   % 设定窗口位置及大小
figure(4);
set(4,'Position',[10,320,300,300]);   % 设定窗口位置及大小
figure(5);
set(5,'Position',[330,320,300,300]);  % 设定窗口位置及大小
figure(6);
set(6,'Position',[650,320,300,300]);  % 设定窗口位置及大小
figure(7);
set(7,'Position',[35,125,300,300]);   % 设定窗口位置及大小
figure(8);
set(8,'Position',[355,125,300,300]);  % 设定窗口位置及大小
figure(9);
set(9,'Position',[675,125,300,300]);   % 设定窗口位置及大小
figure(10);
set(10,'Position',[35,345,300,300]);   % 设定窗口位置及大小
figure(11);
set(11,'Position',[355,345,300,300]);  % 设定窗口位置及大小
figure(12);
set(12,'Position',[675,345,300,300]);  % 设定窗口位置及大小
figure(13);
set(13,'Position',[60,150,300,300]);   % 设定窗口位置及大小
figure(14);
set(14,'Position',[380,150,300,300]);  % 设定窗口位置及大小
figure(15);
set(15,'Position',[700,150,300,300]);   % 设定窗口位置及大小
figure(16);
set(16,'Position',[60,370,300,300]);   % 设定窗口位置及大小
%figure(11);
%set(11,'Position',[380,370,300,300]);  % 设定窗口位置及大小
%figure(12);
%set(12,'Position',[700,370,300,300]);  % 设定窗口位置及大小

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
figure(1)
plot(t,s);
grid on             % 打开网格
axis([-5*Tb,5*Tb,-1.2,1.2])   % 限定范围
title('原始I双极性不归零码')
xlabel('t (us)')
ylabel('x(t)  (V)')

%原始双极性不归零码的功率谱
P=t2f(s);
P=P.*conj(P)/T;
EP=(EP*(j-1)+P)/j;       %RZ功率谱的累计平均
end
Ps=10*log10(EP+eps);% 换算为db值
figure(2)
plot(f,Ps)
grid on             % 打开网格
axis([-5*Rb,5*Rb,-50,10])   % 限定范围
title('原始双极性不归零码的功率谱')
xlabel('f (MHz)')
ylabel('H(f)  (dB/MHz)')

%高斯滤波器的传递函数
Bb=Tb/0.3;
alpha=sqrt(logm(2)/2/Bb^2);
H=exp(-alpha^2*f.^2);
figure(3)
plot(f,H);
grid on             % 打开网格
axis([-10,10,-0.2,1.2])   % 限定范围
title('高斯滤波器的传递函数')
xlabel('f (MHz)')
ylabel('H(f)  (dB/MHz)')

%预编码
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
figure(4)
plot(t,sa);
grid on             % 打开网格
axis([-5*Tb,5*Tb,-1.2,1.2])   % 限定范围
title('预编码后的双极性不归零码')
xlabel('t (us)')
ylabel('sa(t)  (V)')

%高斯滤波后的双极性不归零码
figure(5)
send=real(f2t(t2f(s).*H));
plot(t,send)
grid on             % 打开网格
axis([-5*Tb,5*Tb,-1.2,1.2])   % 限定范围
title('高斯滤波后的双极性不归零码')
xlabel('t (us)')
ylabel('y(t)  (V)')

%串并转换

%It
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
figure(6)
plot(t,It);
grid on;
axis([-5*Tb,5*Tb,-1.2,1.2])   % 限定范围
for k=N:-1:L+1,It(k)=It(k-L);end
figure(7)
plot(t,It);
grid on;
axis([-5*Tb,5*Tb,-1.2,1.2])   % 限定范围
title('It')
xlabel('t (us)')
ylabel('It (V)')

%Qt
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
figure(8)
plot(t,Qt);
grid on;
axis([-5*Tb,5*Tb,-1.2,1.2])   % 限定范围
title('Qt')
xlabel('t (us)')
ylabel('Qt (V)')

%Itt
Itt=It.*cos(pi*t/2/Tb);
figure(9)
plot(t,Itt);
grid on             % 打开网格
axis([-5*Tb,5*Tb,-1.2,1.2])   % 限定范围
title('Itt')
xlabel('t (us)')
ylabel('Itt (V)')

%Qtt
Qtt=Qt.*sin(pi*t/2/Tb);
figure(10)
plot(t,Qtt);
grid on             % 打开网格
axis([-5*Tb,5*Tb,-1.2,1.2])   % 限定范围
title('Qtt')
xlabel('t (us)')
ylabel('Qtt (V)')

%GMSK时域波形
gmsk=Itt.*cos(2*pi*fc*t)-Qtt.*sin(2*pi*fc*t);
figure(11)
plot(t,gmsk);
grid on             % 打开网格
axis([-5*Tb,5*Tb,-1.2,1.2])   % 限定范围
title('GMSK时域信号')
xlabel('t (us)')
ylabel('gmsk (V)')

%GMSK功率谱
PP=t2f(gmsk);
Pa=PP.*conj(PP)/T;
EPg=(EPg*(ii-1)+Pa)/ii;       %RZ功率谱的累计平均
end
Pgmsk=10*log10(EPg+eps);% 换算为db值
figure(12)
plot(f,Pgmsk)
grid on             % 打开网格
axis([-2*fc,2*fc,-70,0])   % 限定范围
title('GMSK功率谱')
xlabel('f (MHz)')
ylabel('H(f) (dB/MHz)')

%接收端
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
figure(13)
plot(f,LPF)
grid on             % 打开网格
axis([-2*fc,2*fc,-0.5,1.5])   % 限定范围
title('接收端的低通滤波器')
xlabel('f (MHz)')

%接收端上支路LPF的输出，与Itt相似（图9）
RI=r.*cos(2*pi*fc*t);
RI=real(f2t(t2f(RI).*LPF));
figure(14)
plot(t,RI);
grid on             % 打开网格
axis([-5*Tb,5*Tb,-1.2,1.2])   % 限定范围
title('接收端上支路LPF的输出')
xlabel('t (us)')
ylabel('x(t)  (V)')

%接收端下支路LPF的输出,与Qtt相似（图10）
RQ=-r.*sin(2*pi*fc*t);
RQ=real(f2t(t2f(RQ).*LPF));
figure(15)
plot(t,RQ);
grid on             % 打开网格
axis([-5*Tb,5*Tb,-1.2,1.2])   % 限定范围
title('接收端下支路LPF的输出')
xlabel('t (us)')
ylabel('x(t)  (V)')

%取样
RIt=RI(2*L:2*L:N);   % 上支路取样,每2Tb取样一个值，在整数处取样,-31-32,共M/2=64个
RQt=RQ(L:2*L:N);   % 下支路取样,每2Tb取样一个值，在整数+0.5处取样,-31.5-31.5,共M/2=64个

%码型串并转换
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
figure(16)
plot(t,sy);
grid on             % 打开网格
axis([-5*Tb,5*Tb,-1.2,1.2])   % 限定范围
title('判决后的双极性不归零码')
xlabel('t (us)')
ylabel('sy(t)  (V)')