%###### 基于DVB-T標準,COFDM調制系統的信道估計的原始模型 ######
%######  仿真條件：4個OFDM符號，QPSK調制，二維內插，信道模型  ######
% IFFT_bin_length: IFFT和FFT的點數
% carrier_count: 子載波個數
% bits_per_symbol: 每個載波含有的比特數
% symbols_per_carrier: 不重復的OFDM符號數
% Ｘ：欲發送的二進制比特流

clear all;
close all;

IFFT_bin_length=2048;
carrier_count=1705;
bits_per_symbol=2;
symbols_per_carrier=4;
LI=12;                  % 導頻之間的間隔
GI=2048/4;              % 保護間隔長度
N_number=carrier_count*symbols_per_carrier*bits_per_symbol;

%--------------------------------------------------------------------------
%矩陣初始化
X=zeros(1,N_number);
X1=[];
X2=[];
X3=[];
X4=[];
X5=[];
X6=[];
X7=[];
Y1=[];
Y2=[];
Y3=[];
Y4=[];
Y5=[];
Y6=[];
Y7=[];
XX=zeros(1,N_number);
dif_bit=zeros(1,N_number);
dif_bit1=zeros(1,N_number);
dif_bit2=zeros(1,N_number);
dif_bit3=zeros(1,N_number);
X=randint(1,N_number);           %產生二進制隨機序列（非0即1）

%--------------------------------------------------------------------------
%% QPSK調制：(1 1)->pi/4;(0 1)->3*pi/4;(0 0)->-3*pi/4;(1,0)->-pi/4; %%
s=1-X.*2;
sreal=s(1:2:N_number);
simage=s(2:2:N_number);
X1=sreal+j.*simage;

%--------------------------------------------------------------------------
%% 串并轉換 %%
X2=reshape(X1,carrier_count,symbols_per_carrier).';

%--------------------------------------------------------------------------
%產生隨機導頻信號
Kmin=0;
Kmax=1704;
pilot(:,1)=Kmin;               %使第一列是導頻
for l=0:symbols_per_carrier-1
    p=0;
    k=0;
    while k<=Kmax-12
       k=Kmin+3*mod(l,4)+12*p;
       pilot(l+1,p+2)=k;       %分散導頻位置
       p=p+1;
    end
end
pilot(:,p+2)=Kmax;             %使最后一列也是導頻
pilot=pilot+1;

g=bin2dec('100000000101');     %偽隨機二進制序列生成多項式:x11+x2+1
state=bin2dec('11111111111');  %偽隨機二進制序列生成寄存器初始狀態
N=2^11-1;                      %生成二進制序列長度
train_sym=zeros(symbols_per_carrier,carrier_count);
m=mgen(g,state,N);             %調用mgen函數，生成二進制隨機序列
for i=1:l+1
    train_sym(i,:)=m(1:carrier_count);
end
for i=1:l+1
    train_sym(i,pilot(i,:))=(4/3)*2.*(1/2-train_sym(i,pilot(i,:))); %導頻值
end

%--------------------------------------------------------------------------
%% 插入導頻 %%
signal=1:carrier_count;
X3(:,signal)=X2;
for i=1:l+1
    X3(i,pilot(i,:))=train_sym(i,pilot(i,:));
end

%--------------------------------------------------------------------------
%% IFFT變換 %%
IFFT_modulation=zeros(symbols_per_carrier,IFFT_bin_length);
IFFT_modulation(:,signal)=X3;
X4=ifft(IFFT_modulation,IFFT_bin_length,2);

%--------------------------------------------------------------------------
%% 加循環前綴保護間隔 %%
for k=1:symbols_per_carrier
   for i=1:IFFT_bin_length
      X6(k,i+GI)=X4(k,i);
   end
   for i=1:GI
      X6(k,i)=X4(k,i+IFFT_bin_length-GI);    
   end
end

%--------------------------------------------------------------------------
%% 并串轉換 %%
X7=reshape(X6.',1,symbols_per_carrier*(IFFT_bin_length+GI));

%--------------------------------------------------------------------------
%信道模型:帶多普勒頻移的瑞利衰落信道
fd=200;                           %多普勒頻移
r=6;                              %多徑數
a=[0.123 0.3 0.4 0.5 0.7 0.8];    %多徑的幅度
d=[20 30 40 50 90 130];           %各徑的延遲
T=1;                              %系統采樣周期
th=[90 0 72 144 216 288]*pi./180; %相移
h=zeros(1,carrier_count);
hh=[];
for k=1:r
    h1=a(k)*exp(j*((2*pi*T*fd*d(k)/carrier_count)));
    hh=[hh,h1];
end
h(d+1)=hh;                        %信道沖激響應
H0=fft(h);

channel1=zeros(size(X7));
channel1(1+d(1):length(X7))=hh(1)*X7(1:length(X7)-d(1));
channel2=zeros(size(X7));
channel2(1+d(2):length(X7))=hh(2)*X7(1:length(X7)-d(2));
channel3=zeros(size(X7));
channel3(1+d(3):length(X7))=hh(3)*X7(1:length(X7)-d(3));
channel4=zeros(size(X7));
channel4(1+d(4):length(X7))=hh(4)*X7(1:length(X7)-d(4));
channel5=zeros(size(X7));
channel5(1+d(5):length(X7))=hh(5)*X7(1:length(X7)-d(5));
channel6=zeros(size(X7));
channel6(1+d(6):length(X7))=hh(6)*X7(1:length(X7)-d(6));

Tx_data=X7+channel1+channel2+channel3+channel4;

%--------------------------------------------------------------------------
%% 加高斯白噪聲 %%  
    snr_db=15;                                          %信噪比
    code_power=[norm(Tx_data)]^2/(length(Tx_data));     %信號的符號功率
    bit_power=code_power/bits_per_symbol;               %比特功率 
    noise_power=10*log10((bit_power/(10^(snr_db/10)))); %噪聲功率
    noise=wgn(1,length(Tx_data),noise_power,'complex'); %產生GAUSS白噪聲信號
    Y7=Tx_data+noise;

%--------------------------------------------------------------------------
%% 串并變換 %%
    Y6=reshape(Y7,IFFT_bin_length+GI,symbols_per_carrier).';
    
%--------------------------------------------------------------------------
%% 去保護間隔 %%
    for k=1:symbols_per_carrier
       for i=1:IFFT_bin_length
           Y5(k,i)=Y6(k,i+GI);
       end
    end

%--------------------------------------------------------------------------
%% FFT變換 %%
    Y4=fft(Y5,IFFT_bin_length,2);
    Y3=Y4(:,signal);
  
%--------------------------------------------------------------------------
%% 信道估計 %%
    H=zeros(size(Y3));
    for i=1:l+1 
        H(i,pilot(i,:))=Y3(i,pilot(i,:))./X3(i,pilot(i,:));  %導頻處信道響應
    end
    
%--------------------------------------------------------------------------
%% 內插 %%
    %導頻間隔為12的位置處的內插
    for i=1:l+1     
        for j=pilot(i,2):12:(pilot(i,p+2)-12)
            for n=1:11
                H(i,j+n)=(1-n/12)*H(i,j)+n/12*H(i,j+12);
            end
        end
    end
    %導頻間隔不為12的位置的內插
    i=2;
    for n=1:2
        H(i,pilot(i,1)+n)=(1-n/3)*H(i,j)+n/3*H(i,pilot(i,2));
    end
    for n=1:8
        H(i,pilot(i,p+1)+n)=(1-n/9)*H(i,pilot(i,p+1))+n/9*H(i,pilot(i,p+2));
    end
    
    i=3;
    for n=1:5
        H(i,pilot(i,1)+n)=(1-n/6)*H(i,j)+n/6*H(i,pilot(i,2));
    end
    for n=1:5
        H(i,pilot(i,p+1)+n)=(1-n/6)*H(i,pilot(i,p+1))+n/6*H(i,pilot(i,p+2));
    end
    
    i=4;
    for n=1:8
        H(i,pilot(i,1)+n)=(1-n/9)*H(i,j)+n/9*H(i,pilot(i,2));
    end
    for n=1:2
        H(i,pilot(i,p+1)+n)=(1-n/3)*H(i,pilot(i,p+1))+n/3*H(i,pilot(i,p+2));
    end
    
%--------------------------------------------------------------------------
%% 濾波 %%
    IFFT_H=ifft(H,carrier_count,2);
    h_filter=[IFFT_H(:,1:GI),zeros(symbols_per_carrier,carrier_count-GI)];
    H_ce=fft(h_filter,carrier_count,2);
    
%--------------------------------------------------------------------------
%% 均衡 %%
    Y2=Y3./H_ce;
             
%--------------------------------------------------------------------------
%% 并串變換 %%
    Y1=reshape(Y2.',1,carrier_count*symbols_per_carrier);
    Y0=reshape(Y3.',1,carrier_count*symbols_per_carrier);
    H1=reshape(H_ce.',1,carrier_count*symbols_per_carrier);
    H1=H1(1:carrier_count);%信道響應
    
%--------------------------------------------------------------------------
%% 作圖 %%
    figure;
    plot(abs(H0),'r:');
    hold on;
    plot(abs(H1));
    hold on;
    legend('No Channel Estimation','Channel Estimation');
    title('信道估計前后傳輸函數比較');
    hold off;	
   
%--------------------------------------------------------------------------
%% QPSK解調 %%
    y_re1=sign(real(Y1));
    y_im1=sign(imag(Y1));
    
    y_re0=sign(real(Y0));
    y_im0=sign(imag(Y0));
 
    r1=zeros(1,N_number);
    r1(1:2:N_number)=y_re1;
    r1(2:2:N_number)=y_im1;
    
    r0=zeros(1,N_number);
    r0(1:2:N_number)=y_re0;
    r0(2:2:N_number)=y_im0;
    
    Y=(1-r01.*sqrt(2))./2;            %輸出二進制序列

% %%%%% 導頻位置及導頻個數統計程序段 %%%%%
% % 得到每個OFDM符號的分散導頻位置
% Kmin=0;
% Kmax=1704;
% for l=0:63
%     p=0;
%     k=0;
%     while k<=Kmax-12
%        k=Kmin+3*mod(l,4)+12*p;
%        m(l+1,p+1)=k;
%        p=p+1;
%     end
% end
% 
% % 比較分散導頻位置和連續導頻位置，計算重復位置個數
% n=0;
% pilot=[0 48 54 87 141 156 192 201 255 279 282 333 432 450 483 525 531 618 636 714 759 765 780 804 873 888 918 939 942 969 984 1050 1101 1107 1110 1137 1140 1146 1206 1269 1323 1377 1491 1683 1704];
% for j=1:45
%     pil=pilot(j);
%     for i=1:p
%         if m(2,i)==pil;
%             n=n+1;
%         end
%     end
% end