%繪制調制波形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);                     %仿真時間值
Ak=[1 0 1 0 1 0 1 0];               %產生8個基帶信號
Ak=2*Ak-1;
gt=ones(1,B_sample);         %每碼元對應的載波信號
Akk=sigexpand(Ak,B_sample);   %碼元擴展
temp=conv(Akk,gt);                 %碼元擴展
Akk=temp(1:length(Akk));            %碼元擴展

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
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);
subplot(311)
plot(t/Tb,Akk);
axis([0 8 -1.5 1.5]);
title('基帶波形');

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

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

%解調
for n=1:512;
    if n<=B_sample
        Alpha1(n)=0;
    else Alpha1(n)=Alpha(n-B_sample);
    end
end
a=[0 1 1 1 1 1 1 1 ]
ak=sigexpand(a,B_sample);   %碼元擴展
temp=conv(ak,gt);                 %碼元擴展
ak=temp(1:length(ak));
S_Gmsk1=cos(2*pi*Fc*(t-Tb)+Alpha1+pi/2).*ak; %延遲1bt,移相pi/2
figure
subplot(311)
plot(t/Tb,S_Gmsk1);
axis([0 8 -1.5 1.5]);
title('延遲1bt,移相pi/2GMSK波形');

xt=S_Gmsk1.*S_Gmsk;
x=0;
subplot(312)
plot(t/Tb,xt,t/Tb,x,'r:');
axis([0 8 -1.5 1.5]);
title('相乘后波形');
%低通濾波
Fs=10000;
rp=3;rs=50;
wp=2*pi*50;ws=2*pi*800;
[n,wn]=buttord(wp,ws,rp,rs,'s')
[z,p,k]=buttap(n);  %把濾波器零極點模型轉化為傳遞函數模型
[bs,as]=lp2lp(bp,ap,wn);  %把模擬濾波器的模型轉化為截止頻率為WN的低通濾波器
[b,a]=bilinear(bs,as,Fs)    %用雙線性變換法實現從模到數的轉化
y=filter(b,a,xt);
subplot(313)
plot(t/Tb,y,t/Tb,x,'r:');
axis([0 8 -1.5 1.5]);
title('經過低通濾波器后波形');

for i=1:8
    if y(i*B_sample)>0
        bt(i)=1
    else
        bt(i)=0
    end
end
bt=2*bt-1;
btt=sigexpand(bt,B_sample);   %碼元擴展
temp1=conv(btt,gt);                 %碼元擴展
btt=temp1(1:length(btt));            %碼元擴展
figure
subplot(311)
plot(bt)
title('抽樣值');
axis([0 8 -1.5 1.5]);
subplot(312)
plot(t/Tb,Akk);
axis([0 8 -1.5 1.5]);
title('原基帶波形');
subplot(313)
plot(t/Tb,btt);
axis([0 8 -1.5 1.5]);
title('解調后波形');