clc;clear;
scmpar=scmparset;
linkpar=linkparset(1); % 10 links
%linkpar.MsVelocity=0;
scmpar.NumMsElements=1;
scmpar.NumBsElements=1;
antpar=antparset;
[H delays out]=scm(scmpar,linkpar,antpar);

% Initialize the parameters

NumLoop = 120;%幀長
NumSubc = 128;%IFFT和FFT的點數
carrier_count=100;%子載波個數
NumCP =8;%循環前綴
TI=7 ; %導頻之間的間隔,梳狀導頻
Np=ceil(NumLoop/TI)+1;%導頻數 %加1的原因：使最后一列也是導頻
carriers=1:carrier_count;

 delta_t = out.delta_t;
 f_sym = 3.84e6;                                       % sample frequency
 osr = floor((1/f_sym)/scmpar.DelaySamplingInterval);  % over sample rate
 ts = 1/(f_sym*osr);                                   % sample time_step
 h_delay = round(delays./ts)+1;                        % multipath delay array: samples
 u = scmpar.NumMsElements;
 s = scmpar.NumBsElements;
 n = scmpar.NumPaths;
 h_tmp = zeros(1,NumSubc*osr);
 h = zeros(s,u,NumSubc);
 for ss = 1:s
     for uu = 1:u
             h_tmp(h_delay(1:n)) = H(uu,ss,1:n,1);
             h(ss,uu,:) = resample(h_tmp,1,osr); % downsample by osr to
     end
 end
 h11 = squeeze(h(1,1,:)); %H11 = fft(h11);
 %h21 = squeeze(h(1,2,:)); H21 = fft(h21).';
 %h12 = squeeze(h(2,1,:)); H12 = fft(h12);
 %h22 = squeeze(h(2,2,:)); H22 = fft(h22).';



% Generate the random binary stream for transmit test
BitsTx = floor(rand(1,2*NumLoop*carrier_count)*2);

% Modulate_signal (generates QAM symbols)
Signal0 = reshape(BitsTx,NumLoop*carrier_count,2);  
Signal1 = bi2de(Signal0,2,'left-msb');        %把二進制數轉換成十進制數
QAMTable = [-1-i -1+i 1-i 1+i];
X1 = QAMTable(Signal1+1)/2;   

% Modulate_training (generates QAM symbols)
training_symbols0 = floor(rand(1,2*carrier_count)*2);
training_symbols1 = reshape(training_symbols0,carrier_count,2);
training_symbols2 = bi2de(training_symbols1,2,'left-msb');
training_symbols3 = QAMTable(training_symbols2+1)/2;
training_symbols3 = training_symbols3';
training_symbols=repmat(training_symbols3,1,Np);

%disp(training_symbols)
pilot=1:TI+1:NumLoop+Np;
if length(pilot)~=Np 
    pilot=[pilot,NumLoop+Np];
end
%--------------------------------------------------------
%串并轉換
X2=reshape(X1,carrier_count,NumLoop);
%---------------------------------------------------------
%插入導頻
signal=1:NumLoop+Np;
signal(pilot)=[];
X3(:,pilot)=training_symbols;
X3(:,signal)=X2;

IFFT_modulation=zeros(NumSubc,NumLoop+Np);
IFFT_modulation(carriers,:)=X3;
%IFFT_modulation(:,conjugate_carriers)=conj(X3);
X4=ifft(IFFT_modulation,NumSubc,1);

% Add cyclic prefix
NumAddPrefix = NumSubc+ NumCP;
RowPrefix = (NumSubc - NumCP + 1):NumSubc;  % 復制
X5 = [X4(RowPrefix,:);X4];

%并串轉換
X6=X5(:);

X7 = filter(h11,1,X6);

snr=0:2:40;
M=length(snr);
ber_n=zeros(M,1);
ber=zeros(M,1);
for i=1:M 
    
    code_power=0;
    code_power=[norm(X7)]^2/(length(X7));%信號的符號功率
    bit_power=code_power/2;%比特功率 
    noise_power=10*log10((bit_power/(10^(snr(i)/10))));%噪聲功率(db)
    noise=wgn(length(X7),1,noise_power,'complex');%產生GAUSS白噪聲信號
    
    Y1=X7+noise;
    
%串并變換,去保護間隔
    Y2 = reshape(Y1,NumAddPrefix,NumLoop+Np);
    Y3 = Y2((NumCP+1):NumAddPrefix,:);
    Y4 = fft(Y3,NumSubc,1);
    Y5=Y4(carriers,:);
   
%LS信道估計
  H=[];
  Y6=Y5(:,signal);
  Rx_training_symbols=Y5(:,pilot);
  Rx_training_symbols0=reshape(Rx_training_symbols,carrier_count*Np,1);
  
  training_symbol0=reshape(training_symbols,1,carrier_count*Np);
  training_symbol1=diag(training_symbol0);
  %disp(training_symbols)
  training_symbol2=inv(training_symbol1);
  Hls=training_symbol2*Rx_training_symbols0;  
  Hls1=reshape(Hls,carrier_count,Np);
  HLs=[];
  HLs1=[];
 if ceil(NumLoop/TI)==NumLoop/TI
     for k=1:Np-1
        HLs2=[];
           for t=1:TI
           HLs1(:,1)=(Hls1(:,k+1)-Hls1(:,k))*(t-1)./LI+Hls1(:,k);
           
           HLs2=[HLs2 HLs1];
           end
       HLs=[HLs HLs2];
    end
else
    for k=1:Np-2
        HLs2=[];
        for t=1:TI
           HLs1(:,1)=(Hls1(:,k+1)-Hls1(:,k))*(t-1)./TI+Hls1(:,k);
           HLs2=[HLs2 HLs1];
        end
       HLs=[HLs HLs2];
    end
    HLs3=[];
    for t=1:mod(NumLoop,TI)
        HLs1(:,1)=(Hls1(:,Np)-Hls1(:,Np-1))*(t-1)./TI+Hls1(:,Np-1);
        HLs3=[HLs3 HLs1];
    end;
    HLs=[HLs HLs3];
end
  %Hls1=Hls.';
  %H=repmat(Hls1,symbols_per_carrier,1);%將導頻擴展成symbols_per_carrier*carrier_count矩陣
  R_Data=Y6./HLs;
  Y7=Y6(:);
  Y8 = zeros(NumLoop*carrier_count,2);
  for m = 1:NumLoop*carrier_count
 
    Real = real(Y7(m));
    Imag = imag(Y7(m));
    if( abs((Real -1)) < abs((Real +1  )) )
        Y8(m,1) = 1;
    else
         Y8(m,1) = 0;
     end
     if( abs((Imag -1)) < abs((Imag +1  )) )
         Y8(m,2) = 1;
     else
         Y8(m,2) = 0;
     end
  end
Y9=Y8(:)';
R_Data=R_Data(:);
SymDec = zeros(NumLoop*carrier_count,2);

for m = 1:NumLoop*carrier_count
 
    Real = real(R_Data(m));
    Imag = imag(R_Data(m));
    if( abs((Real -1)) < abs((Real +1  )) )
        SymDec(m,1) = 1;
    else
         SymDec(m,1) = 0;
     end
     if( abs((Imag -1)) < abs((Imag +1  )) )
         SymDec(m,2) = 1;
     else
         SymDec(m,2) = 0;
     end
end
BitsRx=SymDec(:)';

% Bit Error
[Num_n,Ber_n] = symerr(BitsTx,Y9);
ber_n(i)=Ber_n;
[Num,Ber] = symerr(BitsTx,BitsRx);
ber(i)=Ber;
end
figure(1)
axis([0,40,0,1]);
semilogy(snr,ber_n,'b');
hold on
semilogy(snr,ber,'r');
xlabel('SNR(dB)');
ylabel('BER');
title('ls vs no estimation');
legend('none','ls');
hold on
grid on
ber_n
ber