%圖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)')