%GMSK調制，維特比解調()
clear all
close all

%--------------------------------------------------------------------------
L = 3;                                  %關聯長度
state_number = 2^L*4;                   %計算狀態數
phase_state = [0 pi/2 pi 3*pi/2];       %所有相位狀態
%**************************************************************************

%--------------------------------------------------------------------------
%GMSK信號的所有狀態
state_all = state_all_build(L,phase_state);
%**************************************************************************

%--------------------------------------------------------------------------
%GMSK信號的狀態轉移矩陣
state_trans_matrix = state_trans_build(state_number,L,state_all);
%**************************************************************************

%--------------------------------------------------------------------------
%統計轉入到每個狀態中的狀態序號
state_in = state_in_build(state_number,state_trans_matrix);
%**************************************************************************

%--------------------------------------------------------------------------
%參數設置
data_len = 30000;          %碼元個數
sample_number = 8;          %采樣個數
Rb = 24000;                 %碼元速率
fc = 96000;                 %載波頻率
alpha = 0.25;               %BbTb值
%**************************************************************************

%--------------------------------------------------------------------------
%隨機產生傳輸信號
data = rand_binary(data_len);
%**************************************************************************

%--------------------------------------------------------------------------
%GMSK基帶調制
[signal_out,I_out,Q_out,phase,gt,qt] = mod_gmsk(data,data_len,sample_number,Rb,alpha);
%**************************************************************************

%--------------------------------------------------------------------------
%中頻搬移
I_temp=interp(I_out,4);
Q_temp=interp(Q_out,4);

fc=96000;
Fs=fc*sample_number;
t=1/Fs:1/Fs:length(I_temp)*1/Fs;
signal_i=I_temp.*cos(2*pi*fc*t);
signal_q=Q_temp.*cos(2*pi*fc*t);
signal_mod=I_temp.*cos(2*pi*fc*t)-Q_temp.*sin(2*pi*fc*t);
%**************************************************************************
% [Pxx,w] = pburg(signal_mod(1:1024),34,1024,1*Fs);
% psdplot(Pxx(257:513),w(257:513)-fc,'Hz','','已調CPM信號功率譜密度');           % fc要取在Fs/4(w=[0-Fs/2])
amp=[1 0 0.69 0 0.24];
dj=length(amp);
for i = 1:dj
    signal_mod=signal_mod+amp(i)*delay(signal_mod,i-1,sample_number);
end

%--------------------------------------------------------------------------
%加噪聲
for SNR = 0:15
% signal_mod1 = awgn(signal_mod,snr);

noise=randn(1,length(signal_mod));
b=0.2;c=0.22;d=0.28;e=0.3;
F1=[0,b,c,d,e,1];
A1=[0,0,1,1,0,0];
bpf=firls(300,F1,A1);

noise1=conv(noise,bpf);
noise_out=noise1(150+1:150+length(signal_mod));

Power_noise=sum(noise_out.^2);
Power_signal=sum(signal_mod.^2);
Coe_SNR=sqrt(1*Power_signal/Power_noise/(10^(SNR/10)));

signal_mod1 = signal_mod+(Coe_SNR)*noise_out;    %加入噪聲以后的信號


    %--------------------------------------------------------------------------
    %解調，去載波
    N=300;                                              % 濾波器的階數為(N+1)  
    F=[0,fc-1000,fc+1000,Fs/2]*2/Fs;
    A=[1,1,0,0];
    lpf=firls(N,F,A);
    
    I_dem=signal_mod1.*cos(2*pi*fc*t)*2;
    I_dem=conv(I_dem,lpf);
    I_dem=I_dem(N/2+1:N/2+length(I_temp));
    Q_dem=signal_mod1.*sin(2*pi*fc*t)*2;
    Q_dem=conv(Q_dem,lpf);
    Q_dem=-Q_dem(N/2+1:N/2+length(I_temp));
    
    I_dem_out=zeros(1,length(I_dem)/4);         % 抽取
    Q_dem_out=zeros(1,length(Q_dem)/4);
    for i=1:length(I_dem_out)
       I_dem_out(i)=I_dem(4*(i-1)+1);
       Q_dem_out(i)=Q_dem(4*(i-1)+1);
    end;
    %**************************************************************************

    %--------------------------------------------------------------------------
    %維特比算法
    decision = VA2(I_dem_out,Q_dem_out,data_len,sample_number,state_number,state_all,state_in,L,qt,9);
    %**************************************************************************

    %--------------------------------------------------------------------------
    %計算誤碼率
    [num,ber(SNR+1)] = symerr(data,decision);
    %**************************************************************************
end
%**************************************************************************

%--------------------------------------------------------------------------
%誤碼率曲線
semilogy([0:15],ber,'r*-');
%**************************************************************************