%###### 基于DVB-T標準,COFDM調制系統的利用導頻信號進行符號粗同步 ######
%######  仿真條件：68個OFDM符號，64QAM調制，萊斯信道  ######
% 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=6;
symbols_per_carrier=74; % 總的OFDM符號數(4=10-3-3)
LI=12;                  % 導頻之間的間隔
GI=2048/4;              % 保護間隔長度
N_number=carrier_count*symbols_per_carrier*bits_per_symbol;

%--------------------------------------------------------------------------
%矩陣初始化
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）

%--------------------------------------------------------------------------
%% 64QAM調制 %%
[sreal simage]=QAM_64(X);        %調用64QAM調制函數
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=zeros(symbols_per_carrier,carrier_count);
train_sym=zeros(symbols_per_carrier,carrier_count);
for i=1:l+1
    m=mgen(g,state,N);                      %調用mgen函數，生成二進制隨機序列
    train(i,:)=m(1:carrier_count);
end

for i=1:l+1
    train_sym(i,pilot(i,:))=(4/3)*2.*(1/2-train(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

%% 插入連續導頻 %%
CP=[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]+1;                    %45個連續導頻位置

for i=1:length(CP)
    train_sym(:,CP(i))=(4/3)*2.*(1/2-train(:,CP(i)));%連續導頻值
end
for i=1:length(CP)
    X3(:,CP(i))=train_sym(:,CP(i));                  %插入連續導頻
end

%% 插入傳數參數信令 %%
TPS=[34 50 209 346 413 569 595 688 790 901 1073 1219 1262 1286 1469 1594 1687]+1; %17個傳數參數信令位置

for i=1:length(TPS)
    train_sym(:,TPS(i))=randint;
    train_sym(:,TPS(i))=2*(1/2-train_sym(:,TPS(i))); %傳數參數信令值 
end
for i=1:length(TPS)
    X3(:,TPS(i))=train_sym(:,TPS(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));

% %--------------------------------------------------------------------------
% %% 加入基于DVB-T的萊斯多徑衰落信道(固定接受模式) %%
% [Tx_data h]=dvb_channel(X7);     %調用多徑加入函數
Tx_data=X7;

%--------------------------------------------------------------------------
%% 加高斯白噪聲 %%  
snr_db=15;                                          %信噪比
X2=X2(4:symbols_per_carrier-3,:);                   %取中間連續四個符號
    
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;


N=IFFT_bin_length;
Ng=GI;
r=Y7(20001:30000);           %任意選取一段數據進行估計
% r=Y7;

%--------------------------------------------------------------------------
% 最大似然法
N=IFFT_bin_length;
G=GI;
T=1:3*N+2*G;
for i=1:length(T)            % FFT窗移動的距離
    data1=r(i:G+i-1);        % 取GI長個數據存data1
    data2=r(N+i:N+G+i-1);    % 取相距IFFT_bin_length的GI長個數據存data2
    s(i)=data1*data2';       % 求互相關值    
    cor(i)=abs(s(i)); 
%     cor(i)=abs(real(data1))*(abs(real(data2)))'+abs(imag(data1))*(abs(imag(data2)))';
end

%--------------------------------------------------------------------------
% 作圖
subplot(211);
plot(T,cor);
xlabel('載波數');
ylabel('相關值');
grid on;

[Max Max_i]=max(cor(1:N));                  % Max_i對應一個符號的第一個數據
start_place=Max_i;
for i=1:length(T)
    deltaF(i)=1/(2*pi)*angle(conj(s(i)));
    r1=r(i:G+i-1);
    r2=r(i+N:G+i+N-1);
    e(i)=1/(2*pi)*atan(imag(r1*r2')/(real(r1*r2')));
end
% deltaF=1/(2*pi)*angle(conj(s(Max_i)));      % 分數頻偏估計值
% r1=r(Max_i:G+Max_i-1);
% r2=r(Max_i+N:G+Max_i+N-1);
% e=1/(2*pi)*atan(imag(r1*r2')/(real(r1*r2')));
% disp(start_place);disp(e);
subplot(212);  
plot(T,e);
xlabel('載波數');
ylabel('粗頻偏估計值');
axis(1,length(T),-0.5,0.5);
grid on;