%-------------------------------------------------------------------------%
%                            卷積編碼部分
%-------------------------------------------------------------------------%
%  		cnv_encd(g,k0,input)
%  		determines the output sequence of a binary convolutional encoder
%  		g is the generator matrix of the convolutional code
%        	with n0 rows and l*k0 columns. Its rows are g1,g2,...,gn.
%  		k0 is the number of bits entering the encoder at each clock cycle.
%  		input is the binary input seq.

%  Check to see if extra zero-padding is necessary.
clc
clear all
cnv_g=[1 0 1 1 0 1 1;1 1 1 1 0 0 1;1 1 0 0 1 0 1;1 0 1 1 0 1 1];cnv_q=32*3;
cnv_input=round(rand(1,384*cnv_q));cnv_k0=1;%分別的輸入cnv_q組384位,輸出cnv_q組4*384+24的母碼
if rem(length(cnv_input),cnv_k0) > 0
  cnv_input=[cnv_input,zeros(size(1:k0-rem(length(cnv_input),cnv_k0)))];
end
cnv_n=length(cnv_input)/cnv_k0;
%  Check the size of matrix g.
if rem(size(cnv_g,2),cnv_k0) > 0
  error('Error, g is not of the right size.')
end
%  Determine l and n0.
l=size(cnv_g,2)/cnv_k0;
cnv_n0=size(cnv_g,1);
%  add extra zeros
%  每384位插入l-1個零點
u_tmp=zeros(size(1:(l-1)*cnv_k0));
cnv_u=[u_tmp,cnv_input(1:384)];
for i=1:cnv_q-1
  cnv_u=[cnv_u,u_tmp,cnv_input(384*i+1:384*(i+1))];
end
cnv_u=[cnv_u,u_tmp];
% cnv_u=[u_tmp,cnv_input(1:384),u_tmp,cnv_input(384+1:384*2),u_tmp,cnv_input(384*2+1:384*3),u_tmp,cnv_input(384*3+1:384*4),u_tmp];
%  Generate uu, a matrix whose columns are the contents of 
%  conv. encoder at various clock cycles.
cnv_u1=cnv_u(l*cnv_k0:-1:1);
for i=1:cnv_n+(l-1)*cnv_q-1
  cnv_u1=[cnv_u1,cnv_u((i+l)*cnv_k0:-1:i*cnv_k0+1)];
end
cnv_uu=reshape(cnv_u1,l*cnv_k0,cnv_n+(l-1)*cnv_q);
%  Determine the output
cnv_output=reshape(rem(cnv_g*cnv_uu,2),1,cnv_n0*(cnv_n+(l-1)*cnv_q));
%-------------------------------------------------------------------------%
%                      串并轉換和卷積刪余(puncture)
%刪余原則:采用P=2的保護級別,對于16Kbit/s的數據率,I=384,L=12,L1=5,L2=7,PI1=14,PI2=13,
%刪余輸入端384*4+24.
%-------------------------------------------------------------------------%
VP1=[1 1 1 0 1 1 1 0 1 1 1 0 1 1 0 0 1 1 1 0 1 1 1 0 1 1 1 0 1 1 0 0];
VP2=[1 1 1 0 1 1 1 0 1 1 1 0 1 1 0 0 1 1 1 0 1 1 0 0 1 1 1 0 1 1 0 0];
VT=[1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0];
punc_tmp=(reshape(cnv_output,384*4+24,cnv_q))';
punc_out=zeros(cnv_q,1040);
for i=1:cnv_q
        %將矩陣的每一行分為384*4+24=32*（4*12）+24的子矩陣
    punc_tmptmp=(reshape(punc_tmp(i,1:384*4),32,4*12))';
    punc_VP1=punc_tmptmp(1:20,:);
    punc_VP2=punc_tmptmp(21:48,:);
    punc_VT=punc_tmp(i,384*4+1:384*4+24);
    %若刪余矢量的某元素為零則子矩陣的該列丟棄
    a=0;b=0;c=0;
    for j=1:32
        if VP1(j)==1
            a=a+1;
            punc_VP1tmp1(:,a)=punc_VP1(:,j);
        end
    end
    for j=1:32
        if VP2(j)==1
            b=b+1;
            punc_VP2tmp1(:,b)=punc_VP2(:,j);
        end
    end
    for j=1:24
        if VT(j)==1
            c=c+1;
            punc_VTtmp1(c)=punc_VT(j);
        end
    end
      punc_VP1tmp=(punc_VP1tmp1)';
      punc_VP2tmp=(punc_VP2tmp1)';
      punc_VP1tmp=reshape(punc_VP1tmp,1,20*22);%變行向量
      punc_VP2tmp=reshape(punc_VP2tmp,1,28*21);
      punc_out(i,:)=[punc_VP1tmp punc_VP2tmp punc_VTtmp1];%結果是row 為單位的
end
    punc_outtmp=punc_out';
    punc_outcolumn=punc_outtmp(:);
%-------------------------------------------------------------------------%
%            QPSK映射，頻率交織，差分調制，循環嵌綴，加NULL、PRS,IFFT，信道
%-------------------------------------------------------------------------%
% echo off;clear all;close all;clc;
% fprintf( 'OFDM仿真\n') ;
% tic
% --------------------------------------------- 
%                   參數定義                     
% --------------------------------------------- 
% Initialize the parameters
%test:
%punc_outcolumn=(round(rand(1,1536*8)));
punc_outcolumn=punc_outcolumn';
clear i
NumSubc = 1536;
NumLoop = length(punc_outcolumn)/(NumSubc*2);
NumCP = 384;
SyncDelay = 0;
NumIFFT=2048;
% 子載波數            1536
% 位數/ 符號          2
% 符號數/ 載波        Numloop
% 訓練符號數          
% 循環前綴長度        384   
% 調制方式            QPSK
% 多徑信道數          3 
% IFFT Size          2048  
% 信道最大時延        
% --------------------------------------------- 
%                   QPSK MODULATION              
% --------------------------------------------- 
% Generate the random binary stream for transmit test
BitsTx = punc_outcolumn';
% Modulate (Generates QPSK symbols)
% input: BitsTx(1,NumLoop*NumSubc); output: SymQAM(NumLoop,NumSubc/2)
SymQPSKtmp = reshape(BitsTx,2,NumLoop*NumSubc).';
SymQPSKtmptmp = bi2de(SymQPSKtmp,2,'left-msb');
%--------------------------------------------------------------------
% 函數說明：
% bin2dec(binarystr) interprets the binary string binarystr and returns the
% equivalent decimal number.
%   bi2de是把列向量的每一個元素都由2進制變為10進制
%  D = BI2DE(...,MSBFLAG) uses MSBFLAG to determine the input orientation.
%     MSBFLAG has two possible values, 'right-msb' and 'left-msb'.  Giving a
%     'right-msb' MSBFLAG does not change the function's default behavior.
%   Giving a 'left-msb' MSBFLAG flips the input orientation such that the
%     MSB is on the left.
% % %   D = BI2DE(...,P) converts a base P vector to a decimal value.
% %     Examples:
% %     >> B = [0 0 1 1; 1 0 1 0];
% %     >> T = [0 1 1; 2 1 0];
% %    >> D = bi2de(B)     >> D = bi2de(B,'left-msb')     >> D = bi2de(T,3)
% %     D =                 D =                            D =
% %         12                   3                             12
% %          5                  10                              5
%--------------------------------------------------------------------
% QPSK modulation 
% 00->1+i,01->1-i,10->-1+i,11->-1-i
% 利用查表法進行QPSK星座映射
QPSKTable = [1+i 1-i -1+i -1-i];
SymQPSK = QPSKTable(SymQPSKtmptmp+1);
%--------------------------------------------------------------------------
%                                頻率交織
%--------------------------------------------------------------------------
SymIFFT=reshape(SymQPSK,NumSubc,NumLoop);
A=zeros(1,NumIFFT);
 for ii=2:NumIFFT
     A(ii)=mod(13*A(ii-1)+511,NumIFFT);
 end
 D=zeros(1,NumSubc);
 ii=1;
 for kk=1:NumIFFT
     if A(kk)<=1792&&A(kk)>=256&&A(kk)~=1024
         D(ii)=A(kk);
         ii=ii+1;    
     end
 end
 interleavdata=zeros(NumSubc+1,NumLoop);
 D1=D-255;
 for kk=1:NumSubc
     interleavdata(D1(kk),:)=SymIFFT(kk,:);  
 end
 %-------------------------------------------------------------------------
 %                              PRS的產生，差分調制
 %-------------------------------------------------------------------------
 h1=[0 2 0 0 0 0 1 1 2 0 0 0 2 2 1 1];h2=[0 3 2 3 0 1 3 0 2 1 2 3 2 3 3 0];
 h3=[0 0 0 2 0 2 1 3 2 2 0 2 2 0 1 3];h4=[0 1 2 1 0 3 3 2 2 3 2 1 2 1 3 2];
 Htmp=[h1;h2;h3;h4];
 H=[Htmp Htmp];
 ii=1:4;
 n1=[1 2 0 1 3 2 2 3 2 1 2 3 1 2 3 3 2 2 2 1 1 3 1 2];
 n2=[3 1 1 1 2 2 1 0 2 2 3 3 0 2 1 3 3 3 3 0 3 0 1 1];
 N=[n1 n2];
 h=[H;H;H;H;H;H;H;H;H;H;H;H];
 PHItmp=zeros(size(h));
  for kk=1:32
      PHItmp(:,kk)=h(:,kk)+N';
  end
  PHI=pi/2*PHItmp; 
  tmp=PHI';
  PHI1=tmp(:);% 生成phi值
  Z=exp(j*PHI1);
  PRS=Z;
  %生成對應interleavdata的PRS，即在PRS中間插入一位零
  PRS=[PRS(1:NumSubc/2);0;PRS(NumSubc/2+1:NumSubc)];
  %差分調制
  DQpsk=zeros(NumSubc+1,NumLoop);
  DQpsk(:,1)=PRS.*interleavdata(:,1);
  for kk=2:NumLoop
      DQpsk(:,kk)=DQpsk(:,kk-1).*interleavdata(:,kk)
  end
  %加Null和PRS符號
  Null=zeros(NumSubc+1,1);
  DQpsk=[Null PRS DQpsk];
  %解差分調制：y(:,1)=DQpsk(:,1)./(PRS+10^-16);for...
  %------------------------------------------------------------------------
  %                                IFFT
  %------------------------------------------------------------------------
  %插零，組成2048點
  IFFTtmp=[zeros((NumIFFT-NumSubc)/2,NumLoop+2);DQpsk;zeros((NumIFFT-NumSubc)/2-1,NumLoop+2)];
  SymIFFT=ifft(IFFTtmp,NumIFFT);
 % --------------------------------------------- 
 %             Add cyclic prefix                 
 % --------------------------------------------- 
 NumAddPrefix = NumIFFT + NumCP;
 SymCP = zeros(NumAddPrefix,NumLoop+2);
 RowPrefix = (NumIFFT - NumCP + 1):NumIFFT;
 SymCP = [SymIFFT(RowPrefix,:);SymIFFT];
% --------------------------------------------- 
%             Go through the channel            
% --------------------------------------------- 
% %channel
Data_TX=reshape(SymCP,1,(NumIFFT + NumCP)*(NumLoop+2));
d1=15;a1=0.2;
d2=18;a2=0.3;
copy1=a1.*[zeros(1,d1) Data_TX(1:length(Data_TX)-d1)];
copy2=a2.*[zeros(1,d2) Data_TX(1:length(Data_TX)-d2)];
Data_TX=Data_TX+copy1+copy2;%multipath
snr=20;
SymCh = awgn(Data_TX,snr,'measured');% Add the AWGN
% --------------------------------------------- %