%噪聲加在時域上,四種比較：理想信道、導頻、插值，EM
clc;
clear;
echo off;
T = 0.01;             % set the simulation time
Td =1e-6 ;      % set the bit signal rate , etc 10M hz
Ts = Td ;          % set the sample rate
SNRdB=10:20:70;     % the loop of the SNRindB 
Td2Ts = Td/Ts;
dopler = 50;       % set the dopler shift
baseType = 2; % set the baseband mapping Type , 1: BPSK ; 2 : QPSK ; 3 : 16QAM
subcaNUM = 32 ;       % the number of the OFDM subcarrier
fftdm    = subcaNUM ;       % the number of the FFT points
Tx       = 2  ;    % the number of the transfer attena
Rx       = 2 ;     % the number of the receive attena
cp_l     = 1/4;     % the relative length of the cp 
simulNUM = T/Ts ;  % calculate the simulation points
Dist=60;%導頻的間距
Dist_em=Dist/Dist*1;
sig=0.5:0.5:2;
msen=1:158;
Eav = 1;           % QPSK符號的平均能量 = 1+1


bitSource = bitSource_Gen(simulNUM);%產生一系列的01序列


%bitSource = ones(size(bitSource));
dSource = baseMapping(bitSource,baseType); %按照不同的方式進行調制

% ==================the S2P function
P_dSource = P2S( dSource,subcaNUM);

% ==================the STBC encode should be insert here
P_dSource_tr = P_dSource;   % ====the 1TX1Rx;
%[P_dSource_tr1,P_dSource_tr2] = STBC_2Tx1Rx(P_dSource);
%[P_dSource_tr1,P_dSource_tr2]=VblastTrans_2Tx2Rxintp
% ==================the end of STBC encode

% ==================add training symbol here
%[P_dSource_pilot,pilotNUM] = addpilot(P_dSource_tr,Dist);
[P_dSource_pilot,pilotNUM] = addpilot_em(P_dSource_tr,Dist);
size_P_dSource_pilot=size(P_dSource_pilot);%其實是要算P_dSource_pilot的大小
% ==================The end of the adding pilot

ofdm_dSource = ifft(P_dSource_pilot,fftdm);     % ifft comuputation



Td_ofdm = Ts* Td2Ts * subcaNUM;
cp_length = round(fftdm * cp_l);%cp_l保護間隔是加在頻域上的
tempo_Signal_cp = [ofdm_dSource( fftdm - cp_length+1 : fftdm ,:);ofdm_dSource];  % add cp



tmp = size(tempo_Signal_cp);
ofdm_tempo = reshape(tempo_Signal_cp,1,tmp(1)*tmp(2)) + eps;%1e-5;    % 1e-5 is to avoid NaN%重新變成串行的哦
clear tmp;

 % ==============generate the Rayleigh Multipath fading channel
sim('Mph_Rayleigh_channel',Ts*(length(ofdm_tempo)-1));   % ==========the channel is defined in the simulink model "Mph_Rayleigh_channel"
%第二個參數是指時間軸上的長度
Rayleigh_channel = rec_Signal0 ./ofdm_tempo';
% ===============end of the Rayleigh Multipath fading channel

%for n=1:4
%sigma0=sig(3);
for k = 1:length(SNRdB)
    snr_indB =SNRdB(k);
    %if    snr_indB<30 
    %sigma0=1.3;
    %else 
    %sigma0=0.15;
    %end
    
  snr = 10^(snr_indB/10);   
  N0 = Eav/snr;
  if snr_indB<30
      sigma0=sqrt(N0/2);
  elseif snr_indB<50 sigma0 = 2*sqrt (N0/2);    %加入雙邊帶功率譜
  elseif snr_indB<60 sigma0=5*sqrt(N0/2);
  elseif snr_indB<80 sigma0=10*sqrt(N0/2);
  elseif snr_indB<100 sigma0=100*sqrt(N0/2);
  else sigma0=1000*sqrt(N0/2);
  end
    %else sigma0 = 0.015;
    %end
  %if    snr_indB<30 
  %sigma0=1.3;
  %else 
  %sigma0=0.015;
  %end
% =====================noise is added to channel here  
%rec_Signal = awgn(rec_Signal0,snr_indB,'measured');
% ==================then end of the adding noise


% ==================if there is no timing error
%rec_P_tempo = reshape(rec_Signal0,size(tempo_Signal_cp));
rec_P_tempo = reshape(rec_Signal0,size(tempo_Signal_cp));  
%rec_P_tempo = reshape(rec_Signal0,size(tempo_Signal_cp));
rec_tempo_Signal = rec_P_tempo(cp_length +1 : fftdm+cp_length ,:);    % remove cp



rec_freq_Signal_ideal = fft(rec_tempo_Signal,fftdm);         % fft computation
H_ideal=rec_freq_Signal_ideal./P_dSource_pilot;

%rec_freq_Signal=awgn(rec_freq_Signal_ideal,snr_indB,'measured');
rec_freq_Signal=add_noise(rec_freq_Signal_ideal,snr_indB,Eav);
% ==================The start of the channel estimation section

 %if n==1
 %rec_freq_Signal_nopilot = channel_esti_after_em(H_ideal,rec_freq_Signal,Dist,pilotNUM,1);
 %elseif n==2 
 %channel_we_esti = channel_esti_em727m(rec_freq_Signal,Dist,pilotNUM,sigma0,Dist_em);
[channel_we_esti,kh]= channel_esti_em730a(rec_freq_Signal,Dist,pilotNUM,sigma0);
rec_freq_Signal_nopilot = channel_esti_after_em(channel_we_esti,rec_freq_Signal,Dist,pilotNUM,1);
ERR=abs(channel_we_esti-H_ideal);
MSE(k)=sum(sum(ERR))'/(size_P_dSource_pilot(2)*size_P_dSource_pilot(1));
 %elseif n==3
 %rec_freq_Signal_nopilot = channel_esti(rec_freq_Signal,Dist,pilotNUM,1);
 %elseif n==4 
 %rec_freq_Signal_nopilot = channel_esti_intp_time(rec_freq_Signal,Dist,pilotNUM,1);
 %end
% ==================The end of the channel estimation section
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%add something
%H_ideal=rec_freq_Signal_./P_dSource_pilot;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%here

% ===================the STBC decode should be inserted here
  
% ===================end of the STBC decode
rec_freq_Signal_soft = reshape(rec_freq_Signal_nopilot,1,length(rec_freq_Signal_nopilot(1,:)) * length(rec_freq_Signal_nopilot(:,1)) );

%
%F = [1:length(rec_freq_Signal_soft)];
%subplot(2,1,1);
%plot(F(1:32),real(rec_freq_Signal_soft(1:32)));
%subplot(2,1,2);
%plot(F(1:32),imag(rec_freq_Signal_soft(1:32)));


rec_bitSource = ofdmSlice_c(rec_freq_Signal_soft,baseType);
%rec_bitSource = ofdmSlice(rec_freq_Signal_soft,baseType);
BER(k) = bitErrorCalc(bitSource,rec_bitSource);

clc;

     sizekh=size(kh);
     timel=1:sizekh(2);
    % if snr_indB==10
 %    semilogy(timel,kh,'r--x');xlabel('t'),ylabel('iteration number');hold on;
  %    elseif snr_indB==30
  %    semilogy(timel,kh,'b-o');hold on;
   %   elseif snr_indB==50
   %   semilogy(timel,kh,'g-+');hold on;
   %   elseif snr_indB==70
    %  semilogy(timel,kh,'y-*');hold on;
  %   end


    %end                                % =============the end of the SNRindB loop

    %if  n==3 semilogy( SNRdB,BER,'g--o'); xlabel('SNR(dB)'),ylabel('BER');hold on ;
    %elseif n==1 semilogy(SNRdB,BER,'r-*');hold on;
    
    %elseif n==2 semilogy(SNRdB,BER,'b-+');hold on;
    %elseif n==4 semilogy(SNRdB,BER,'y-x');hold on;
         %semilogy( SNRdB,MSE,'r--o'); xlabel('SNR(dB)'),ylabel('MSE');hold on ;
         semilogy( SNRdB,BER,'r--o'); xlabel('SNR(dB)'),ylabel('BER');hold on ; 
         % end
     
 end