% --------------------------------------------- %
% 摘錄自《基于MATLAB的OFDM系統仿真及性能分析》     
% 作者：呂愛琴 田玉敏 朱明華                      %
% 第22卷　第10期  計算機仿真  2005年10月          %
% 整理：  lavabin                              %
% 2006.04.24                          %
% --------------------------------------------- %
% 【問題】 
% 1. 該程序進行了簡單的LMS信道估計，沒有加入與MMSE等其他信道估計算法的比較；
% 2. 程序中循環前綴的設置，作者稱是25％的CP，但是從矩陣的構造并不能看出這一點
%    并且CP為Zero-CP，需要加入與真正復制了數據樣值的循環CP系統進行性能對比；
% 3. 仿真條件為系統處于理想同步情況下。
% 
% 【參考程序】 個人認為該程序參考了Eric Lawrey的學士學位論文程序。
%             原因：二者對每個OFDM符號的載波分配都采用了負頻率對稱排布；
%                  二者對發送數據在載波上的位置分配的設計完全相同；  
%                  調制方式都是DQPSK，調制解調代碼完全一致。

clear all;
close all;
fprintf( '\n OFDM仿真\n \n') ;

% --------------------------------------------- %
%                   參數定義                     %
% --------------------------------------------- %

IFFT_bin_length = 1024
carrier_count   = 200
bits_per_symbol = 2
symbols_per_carrier = 50
% 子載波數            200
% 位數/ 符號          2
% 符號數/ 載波        50
% 訓練符號數          10
% 循環前綴長度        T/4（作者注明，不解為什么是該值）  
% 調制方式            QDPSK
% 多徑信道數          2、3、4（缺省）
% 信道最大時延        7 (單位數據符號)
% 仿真條件            收發之間嚴格同步 

SNR = input('SNR =')
% 輸入信噪比參數

baseband_out_length = carrier_count * symbols_per_carrier * bits_per_symbol
% 計算發送的二進制序列長度

carriers = (1: carrier_count) + (floor(IFFT_bin_length/4) - floor(carrier_count/2))
conjugate_carriers = IFFT_bin_length - carriers + 2
% 構造共軛時間－載波矩陣，以便應用所謂的RCC，Reduced Computational Complexity算法，即ifft之后結果為實數 
% Define the conjugate time-carrier matrix
% 也可以用flipdim函數構造對稱共軛矩陣

% --------------------------------------------- %
%                   信號發射                     %
% --------------------------------------------- %
out = rand(1,baseband_out_length); %隨機數
baseband_out1 = round(out) ; %最接近x的整數
baseband_out2 = floor(out*2) ;%向負無窮取整
baseband_out3 = ceil(out*2)-1 ;%向正無窮取整 
baseband_out4 = randint(1,baseband_out_length);%產生均勻分布的隨機整數矩陣
% 四種生成發送的二進制序列的方法，任取一種產生要發送的二進制序列

if (baseband_out1 == baseband_out2 & baseband_out1 == baseband_out3 )
    fprintf('Transmission Sequence Generated \n \n');
    baseband_out = baseband_out1 ;
else 
    fprintf('Check Code!!!!!!!!!!!!!!!!!!!!! \n \n');
end
% 驗證四種生成發送的二進制序列的方法

convert_matrix = reshape(baseband_out,bits_per_symbol,length(baseband_out)/bits_per_symbol) ;% 改變數組維數、大小
% 函數說明：
% RESHAPE Change size.
%    RESHAPE(X,M,N) returns the M-by-N matrix whose elements
%    are taken columnwise from X.  An error results if X does
%    not have M*N elements.

for k = 1:(length(baseband_out)/bits_per_symbol)
    modulo_baseband(k) = 0 ;
    for i = 1:bits_per_symbol
        modulo_baseband(k) =  modulo_baseband(k) + convert_matrix(i,k)*2^(bits_per_symbol - i) ;
    end
    end
% 每2個比特轉化為整數 0至3
% 采用'left-msb'方式

%-------------------------------------------------------------------------
%  Test by lavabin
%  A built-in function of directly change binary bits into decimal numbers
%-------------------------------------------------------------------------
convert_matrix1 = zeros(length(baseband_out)/bits_per_symbol,bits_per_symbol);%全零數組
convert_matrix1 = convert_matrix' ;
Test_convert_matrix1 = bi2de(convert_matrix1,bits_per_symbol,'left-msb');%二進制轉換為十進制
Test_convert_matrix2 = bi2de(convert_matrix1,bits_per_symbol,'right-msb');
% 函數說明：
%  BI2DE Convert binary vectors to decimal numbers.
%     D = BI2DE(B) converts a binary vector B to a decimal value D. When B
%     is
%     a matrix, the conversion is performed row-wise and the output D is a
%     column vector of decimal values. The default orientation of the binary
%     input is Right-MSB; the first element in B represents the least
%     significant bit.
if (modulo_baseband == Test_convert_matrix1') 
    fprintf('modulo_baseband = Test_convert_matrix1 \n\n\n');
else if (modulo_baseband == Test_convert_matrix2')      
    fprintf('modulo_baseband = Test_convert_matrix2 \n\n\n');
    else
    fprintf('modulo_baseband ~= any Test_convert_matrix \n\n\n'); 
    end
end
% we get the result "modulo_baseband = Test_convert_matrix1".
%-------------------------------------------------------------------------
  
carrier_matrix = reshape(modulo_baseband,carrier_count,symbols_per_carrier)';  
% 生成時間－載波矩陣

% --------------------------------------------- %
%                   QDPSK調制                   %
% --------------------------------------------- %
carrier_matrix = [zeros(1,carrier_count); carrier_matrix];
% 添加一個差分調制的初始相位，為0
for i = 2:(symbols_per_carrier + 1)
    carrier_matrix(i,:) = rem(carrier_matrix(i,:) + carrier_matrix (i-1,:), 2^bits_per_symbol) ;
% 差分調制 rem-求余數
end

carrier_matrix = carrier_matrix*((2*pi)/(2^bits_per_symbol)) ;
% 產生差分相位

[X, Y]=pol2cart(carrier_matrix, ones(size(carrier_matrix,1),size(carrier_matrix,2))); %極或柱坐標變為直角坐標
% 由極坐標向復數坐標轉化 第一參數為相位 第二參數為幅度
% Carrier_matrix contains all the phase information and all the amplitudes
% are the same,‘1’.  
%－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
% 函數說明：
%  POL2CART Transform polar to Cartesian coordinates.
%    [X,Y] = POL2CART(TH,R) transforms corresponding elements of data
%    stored in polar coordinates (angle TH, radius R) to Cartesian
%    coordinates X,Y.  The arrays TH and R must the same size (or
%    either can be scalar).  TH must be in radians.

% [X,Y,Z] = POL2CART(TH,R,Z) transforms corresponding elements of
%    data stored in cylindrical coordinates (angle TH, radius R, height Z) 
%    to Cartesian coordinates X,Y,Z. The arrays TH, R, and Z must be
%    the same size (or any of them can be scalar).  TH must be in radians.
%－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－

complex_carrier_matrix = complex(X, Y) ;% 合成的
%－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
% 函數說明：
%  COMPLEX Construct complex result from real and imaginary parts.
%    C = COMPLEX(A,B) returns the complex result A + Bi, where A and B are
%    identically sized real N-D arrays, matrices, or scalars of the same data type.

% 添加訓練序列
training_symbols = [ 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 ...
-j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 ...
1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 ...
-1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j ...
-1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 ]; 
% 25 times "1 j j 1" 
% 25 times "-1 -j -j -1"
% totally 200 symbols as a row

training_symbols = cat(1, training_symbols, training_symbols) ;%串接成高維數組 
training_symbols = cat(1, training_symbols, training_symbols) ;   
% Production of 4 rows of training_symbols
%-------------------------------------------------------------------------
% 函數說明：
% CAT Concatenate arrays.
%     CAT(DIM,A,B) concatenates the arrays A and B along the dimension DIM.  
%     CAT(2,A,B) is the same as [A,B].
%     CAT(1,A,B) is the same as [A;B].
%------------------------------------------------------------------------- 
%     B = CAT(DIM,A1,A2,A3,A4,...) concatenates the input
%     arrays A1, A2, etc. along the dimension DIM.

complex_carrier_matrix = cat(1, training_symbols, complex_carrier_matrix);
% 訓練序列與數據合并 
% block-type pilot symbols

IFFT_modulation = zeros(4 + symbols_per_carrier + 1,IFFT_bin_length) ;
% % Here a row vector of zeros is between training symbols and data
% symbols!!! 
% % 4 training symbols and 1 zero symbol
% % every OFDM symbol takes a row of "IFFT_modulation" 

IFFT_modulation(: , carriers) = complex_carrier_matrix;
IFFT_modulation(: , conjugate_carriers) = conj(complex_carrier_matrix) ;

%-------------------------------------------------------------------------
%   Test by lavabin  -- Find the indices of zeros 
index_of_zeros = zeros(symbols_per_carrier,IFFT_bin_length - 2*carrier_count);
IFFT_modulation1 = zeros(4 + symbols_per_carrier + 1,IFFT_bin_length);
IFFT_modulation2 = zeros(4 + symbols_per_carrier + 1,IFFT_bin_length);
IFFT_modulation1(6:symbols_per_carrier+5,:) = IFFT_modulation(6:symbols_per_carrier+5,:)==0 ;
for i = 1:symbols_per_carrier
    index_of_zeros(i,:) = find(IFFT_modulation1(i+5,:)==1);%尋找非零元素下標
end
%-------------------------------------------------------------------------

time_wave_matrix = ifft(IFFT_modulation') ; % 進行IFFT操作
time_wave_matrix = time_wave_matrix';
% If X is a matrix, ifft returns the inverse Fourier transform of each column of
% the matrix.

for i = 1: 4 + symbols_per_carrier + 1
    windowed_time_wave_matrix( i, : ) = real(time_wave_matrix( i, : )) ;
end
% get the real part of the result of IFFT
% 這一步可以省略，因為IFFT結果都是實數
% 由此可以看出，只是取了IFFT之后載波上的點，并未進行CP的復制和添加end

ofdm_modulation = reshape(windowed_time_wave_matrix',1, IFFT_bin_length*(4 + symbols_per_carrier + 1) ) ;
% P2S operation
%-------------------------------------------------------------------------