clear all
clc
% 這個程序主要用于仿真兩發一收的STBC-GMSK系統，并用MMSE軟干擾抵消的方法進行SISO迭代檢測
% Copyright Oct.29th 2008, Jun Zhou
% USETC
% for academic use only
      diary  STBC-GMSK_N_2_L_3.txt

% Choose decoding algorithm
dec_alg = input('Please enter the decoding algorithm.(0:Log-MAP: default 0)');
if isempty(dec_alg)
   dec_alg = 0;
end

niter = input('Please enter the maximum number of iteration detection.(default 1)');
if isempty(niter)
    niter = 1;
end

% Frame size
L_total = input('Please enter the frame size(info + tail,default: 512)');
if isempty(L_total)
   L_total = 512;	 % infomation bits plus tail bits
end

% Code generator
g = input('Please enter code generator:(default: g = [1 1 1 1; 1 1 0 1] )');
if isempty(g)
   g = [ 1 1 1 1;
         1 1 0 1];
end

m = input('Please enter the number of modulation order.(power of 2: default 2)');
if isempty(m)
   m = 2;
end

% Select pulse type
pulse = input('Please enter the type of pulse .(0: LRC, 1: GMSK : default 0)');
if isempty(pulse)
   pulse = 0;
end

L = input('Please enter the memory length of the pulse.(L :default 3)');
if isempty(L)
   L = 3;
end

BT = input('Please enter -3dB bandwidth of thd Gaussian pulse.(BT = 0.3, 0.5 :default 0.3)');
if isempty(BT)
   BT = 0.3;
end

% the modulation index
h_m=1;       %molecule of modulation of index
h_n=2;       %numerator of modulation of index
h = h_m/h_n;

Ns = input('Please enter the oversampling rate.(4,8,16: default 16)');
if isempty(Ns)
   Ns = 16;
end

[n,K] = size(g);
md = K - 1;
nstates = 2^md;

% Code rate
rate = 0.5;

% Number of frame errors to count as a stop criterior
berrlim = input('Please enter number of bit errors to terminate:(default 500)');
if isempty(berrlim)
   berrlim = 500;
end

TotalBit = input('Please enter number of simulation bits:(default 100000000)');
if isempty(TotalBit)
   TotalBit = 100000000;
end  

EbN0db = input('Please enter Eb/N0 in dB: (default [5])');
if isempty(EbN0db)
   EbN0db = 5;
end

fprintf('\n\n----------------------------------------------------\n');
if dec_alg == 0
   fprintf(' === Log-MAP decoder === \n');
else
   fprintf(' === Max-Log-MAP decoder === \n'); 
end

if pulse == 0
   fprintf(' === LRC  pulse === \n');
else
   fprintf(' === GMSK pulse === \n'); 
end

fprintf(' Frame size = %6d\n',L_total);
fprintf(' code generator: \n');
for i = 1:n
    for j = 1:K
        fprintf( '%6d', g(i,j));
    end
    fprintf('\n');
end 
fprintf(' The maximum number of iteration = %6d\n',niter);
fprintf(' terminate bit errors = %6d\n', berrlim);
fprintf(' Simulation bit numbers = %6d\n', TotalBit);
fprintf(' The number of modulation order: m = %6d\n', m);
fprintf(' The memory length of the pulse: L = %6d\n', L);
fprintf(' The -3dB bandwidth of thd Gaussian pulse: BT = %f\n', BT);
fprintf(' The modulation index h = %f\n', h);
fprintf(' Eb/N0 (dB) = ');
for i = 1:length(EbN0db)
    fprintf('%10.2f',EbN0db(i));
end
fprintf('\n----------------------------------------------------\n\n');
   
fprintf('+ + + + Please be patient. Wait a while to get the result. + + + +\n');

if pulse == 0 
    q = gen_cpfsk_rc(Ns, L);
else
    q = gen_cpfsk_gmsk(Ns, L, BT);
end
% generate matchaed filter
[mtchd_fltr,D] = gen_mtchd_fltr(Ns, m, L, h, q);
bin_matrix = de2bi([0:m-1],log2(m), 'left-msb');

for nEN = 1:length(EbN0db)
    SNR = EbN0db(nEN) + 10*log10(rate) - 10*log10(Ns) + 10*log10(log2(m));
    snr = 10^(SNR/10);
    sigma = 1/sqrt(2*snr);
    info_sym_num = floor(2*L_total/log2(m));  
    info_bit_num = info_sym_num * log2(m);
    xlast = zeros(1,info_bit_num/2);
    nframe = 0;
    % Clear bit error counter and frame error counter
    nberr(nEN,1:niter) = zeros(1,niter);
    nferr(nEN,1:niter) = zeros(1,niter);
    while lt(nberr(nEN,niter),berrlim) & lt(nframe*(L_total-md),TotalBit)
        nframe = nframe + 1;
        x = round(rand(1, L_total-md));         % info. bits
        en_output = encoderm(x, g) ;            % encoder output (+1/-1)
        [B,ind] = sort(rand(1,info_bit_num));
        mod_input = en_output(ind);             % interleave
        d1 = 2 * mod_input - 1;
        d = reshape(d1,2,info_bit_num/2);
        d_tmp1 = - d(2,:);
        d_tmp2 =   d(1,:);
        d = [d_tmp1;d_tmp2];
        d2 = reshape(d,1,info_bit_num);
        d = [d1;d2];
        alpha = zeros(2,info_bit_num);
        alpha(:,1) = [d1(1);d2(1)];
        for l=1 : 2
            for k = 2:info_bit_num
                alpha(l,k) = (-1)^(k-1) * d(l,k-1) * d(l,k); % Differential precoding
            end
        end
%************************************* Modulation **************************************    
        [src_data1, mod_sig1] = mod_cpfsk_sig(alpha(1,:), 0, h, m, Ns, L, q);
        [src_data2, mod_sig2] = mod_cpfsk_sig(alpha(2,:), 0, h, m, Ns, L, q);
        % Via Rayleigh channel and Add gaussian noise
        N = 2; % 采樣率
        H = sqrt(1/2)*(randn(N,2) + sqrt(-1)*randn(N,2)); % Channel coefficients
        mod_sig = [mod_sig1 mod_sig2].';
        mod_sig = H * mod_sig/sqrt(2);
        noisy_sig = mod_sig + sigma * (randn(size(mod_sig)) + sqrt(-1) * randn(size(mod_sig))); 
        
        %iter_num = niter;
        L_e = zeros(1,info_bit_num);
        for iter = 1:niter
% ************************************* Demodulator ************************************
            L_a = L_e;                     % a priori info.    
            L_all = mmse3(noisy_sig, mtchd_fltr, D, Ns, H, L, sigma, L_a, N, d1); % complete info.
            L_e = (L_all);                 % extrinsic info.
            L_ch(ind) = 0.5 * L_e;         % deinterleave combined channel and extrinsic info.
% *************************************** Decoder **************************************
            L_a = zeros(1,length(L_ch)/2); % a priori info.
            L_all = logmapo(L_ch, g, L_a); % complete info.
            L_e = L_all - 2*L_ch;          % extrinsic info.
            L_e = L_e(ind);                % interleave

% Estimate the info. bits
            xhat = (1-sign(L_all(1:2:info_bit_num)))/2;
% % Number of bit errors in current iteration
%             err(iter) = length(find(xhat(1:L_total-md)~=x));
% Number of bit errors in current iteration
            err(iter) = length(find(xhat(1:L_total-md)~=x));
            if length(find(xhat~=xlast))==0
                err(iter:niter)=err(iter);
                iter_num = iter;
                break;
            end
% Count frame errors for the current iteration
            if gt(err(iter),0)
               nferr(nEN,iter) = nferr(nEN,iter)+1;
            end
            xlast = xhat;
         end % end iter
% Total number of bit errors for all iterations
            nberr(nEN,1:niter) = nberr(nEN,1:niter) + err(1:niter);
            if rem(nframe,100)==0 | nberr(nEN,niter)>=berrlim
               %error_demod/nframe/(L_total-md)
% Bit error rate
               ber(nEN,1:niter) = nberr(nEN,1:niter)/nframe/(L_total-md);
% Frame error rate
               fer(nEN,1:niter) = nferr(nEN,1:niter)/nframe;

% Display intermediate results in process 
               fprintf('************** Eb/N0 = %5.2f db **************\n',EbN0db(nEN));
               fprintf('Frame size = %d rate 1/%d. Modulation order = %d,h = %f\n',L_total,2,m,h);
               fprintf('%d frames transmitted, %d bits in error.\n', nframe, nberr(nEN,niter));
               fprintf('Bit Error Rate (from iteration 1 to iteration %d):\n',niter);
               for i=1:niter
                  fprintf('%8.4e    ', ber(nEN,i));
               end
               fprintf('\n');
               fprintf('Frame Error Rate (from iteration 1 to iteration %d):\n', niter);
               for i=1:niter
                  fprintf('%8.4e    ', fer(nEN,i));
               end
               fprintf('\n');
          
               fprintf('***********************************************\n\n');

% Save intermediate results
                save data EbN0db ber fer 
            end % if 
        
    end % while 
        
end % nEN
 diary off