% 【問題】 
% 1. 該程序進行了簡單的LMS信道估計，沒有加入與MMSE等其他信道估計算法的比較；
% 2. 程序中循環(huán)前綴的設(shè)置，作者稱是25％的CP，但是從矩陣的構(gòu)造并不能看出這一點
%    并且CP為Zero-CP，需要加入與真正復(fù)制了數(shù)據(jù)樣值的循環(huán)CP系統(tǒng)進行性能對比；
% 3. 仿真條件為系統(tǒng)處于理想同步情況下。
% 

clear all;
close all;
fprintf( '\n OFDM仿真\n \n') ;
% --------------------------------------------- %
%                   參數(shù)定義                     %
% --------------------------------------------- %
IFFT_bin_length = 1024;
carrier_count   = 200;
bits_per_symbol = 2;
symbols_per_carrier = 50;
% 子載波數(shù)            200
% 位數(shù)/ 符號          2
% 符號數(shù)/ 載波        50
% 訓(xùn)練符號數(shù)          10
% 循環(huán)前綴長度        T/4（作者注明，不解為什么是該值）  All-zero CP  
% 多徑信道數(shù)          2、3、4（缺省）
% 信道最大時延        7 (單位數(shù)據(jù)符號)
% 仿真條件            收發(fā)之間嚴格同步
SNR =364;
%SNR = input('SNR =') ;% 輸入信噪比參數(shù)
baseband_out_length = carrier_count * symbols_per_carrier * bits_per_symbol;% 計算發(fā)送的二進制序列長度
carriers = (1: carrier_count) + (floor(IFFT_bin_length/4) - floor(carrier_count/2));   % 坐標： (1 to 200) + 156 ,  157 -- 356
conjugate_carriers = IFFT_bin_length - carriers + 2;    % 坐標 ：1024 - (157:356) + 2 = 1026 - (157:356) = (869:670)     
% 構(gòu)造共軛時間－載波矩陣，以便應(yīng)用所謂的RCC，Reduced Computational Complexity算法，即ifft之后結(jié)果為實數(shù) 
% Define the conjugate time-carrier matrix
% 也可以用flipdim函數(shù)構(gòu)造對稱共軛矩陣

% --------------------------------------------- %
%                   信號發(fā)射                     %
% --------------------------------------------- %
out = rand(1,baseband_out_length);
baseband_out = round(out) ;
convert_matrix = reshape(baseband_out,length(baseband_out)/bits_per_symbol,bits_per_symbol) ;
modulo_baseband = bi2de(convert_matrix,bits_per_symbol,'left-msb');%1*10000的十進制矩陣
carrier_matrix = reshape(modulo_baseband,carrier_count,symbols_per_carrier)'; % 生成時間－載波矩陣 50*200

% --------------------------------------------- %
%                   QDPSK調(diào)制                   %
% --------------------------------------------- %
carrier_matrix = [zeros(1,carrier_count); carrier_matrix];% 添加一個差分調(diào)制的初始相位，為0  51*200矩陣
for i = 2:(symbols_per_carrier + 1)
carrier_matrix(i,:) = rem(carrier_matrix(i,:) + carrier_matrix (i-1,:), 2^bits_per_symbol) ;% 差分調(diào)制 
end
carrier_matrix = carrier_matrix*((2*pi)/(2^bits_per_symbol)) ;
% 產(chǎn)生差分相位
[X, Y]=pol2cart(carrier_matrix, ones(size(carrier_matrix,1),size(carrier_matrix,2))); 
% 由極坐標向復(fù)數(shù)坐標轉(zhuǎn)化 第一參數(shù)為相位 第二參數(shù)為幅度
% Carrier_matrix contains all the phase information and all the amplitudes
% are the same,‘1’.  
%－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
% 函數(shù)說明：
%  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) ;%51*200的QPSK調(diào)制信號
figure(1);
plot(complex_carrier_matrix,'*r')
axis([-2,2,-2,2])
grid on
%－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
% 添加訓(xùn)練序列
%－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
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   1*200
training_symbols = [training_symbols;training_symbols] ;  
training_symbols = [training_symbols;training_symbols] ;   %4*200
complex_carrier_matrix = [training_symbols;complex_carrier_matrix] ; % 訓(xùn)練序列與數(shù)據(jù)合并 55*200


IFFT_modulation = zeros(4 + symbols_per_carrier + 1,IFFT_bin_length) ;%55*1024
% % 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;%carriers=[157:356]
IFFT_modulation(: , conjugate_carriers) = conj(complex_carrier_matrix) ;%conjugate_carriers=[869:670]
%-------------------------------------------------------------------------

%-------------------------------------------------------------------------
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.
%IFFT結(jié)果都是實數(shù)
% 由此可以看出，只是取了IFFT之后載波上的點，并未進行CP的復(fù)制和添加end

figure(2);
subplot(4,1,1)
plot(time_wave_matrix)
grid on

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

% --------------------------------------------- %
%                   信道模擬                     %
% --------------------------------------------- %
d1 = 4; a1 = 0.2; d2 = 5; a2 = 0.3; d3 = 6; a3 = 0.4;
d4 = 7; a4 = 0.5;
copy1 = zeros(size(Tx_data)) ;
for i = 1 + d1: length(Tx_data)
copy1(i) = a1*Tx_data( i - d1) ;
end
copy2 = zeros(size(Tx_data) ) ;
for i = 1 + d2: length( Tx_data)
copy2(i) = a2*Tx_data( i - d2) ;
end
copy3 = zeros(size(Tx_data) ) ;
for i = 1 + d3: length(Tx_data)
copy3(i) = a3*Tx_data ( i - d3) ;
end
copy4 = zeros(size(Tx_data) ) ;
for i = 1 + d4: length( Tx_data)
copy4(i) = a4*Tx_data(i - d4) ;
end
Tx_data = Tx_data + copy1 + copy2 + copy3 + copy4;% 4 multi-paths

subplot(4,1,2)
plot(Tx_data)
grid on

Tx_signal_power = var(Tx_data);

%-------------------------------------------------------------------------
% 函數(shù)說明：
%  VAR Variance.
%     For vectors, Y = VAR(X) returns the variance of the values in X.  For
%     matrices, Y is a row vector containing the variance of each column of
%     X.
linear_SNR = 10^( SNR /10) ;
noise_sigma = Tx_signal_power / linear_SNR;
noise_scale_factor = sqrt(noise_sigma) ;
noise = randn(1, length(Tx_data) )*noise_scale_factor;

subplot(4,1,3)
plot(noise)

%-------------------------------------------------------------------------
% 函數(shù)說明：
% Y = randn(m,n) or Y = randn([m n]) returns an m-by-n matrix of random
% entries.
% The randn function generates arrays of random numbers whose elements are
% normally distributed with mean 0 and variance 1.
Rx_Data = Tx_data + noise;
subplot(4,1,4)
plot(Rx_Data)
% --------------------------------------------- %
%                  信號接收                      %
% --------------------------------------------- %
Rx_Data_matrix = reshape(Rx_Data, IFFT_bin_length, 4 + symbols_per_carrier + 1) ;
Rx_spectrum = fft(Rx_Data_matrix) ; 
% %  Suppose precise synchronazition between Tx and Rx
Rx_carriers = Rx_spectrum( carriers, : )';
Rx_training_symbols = Rx_carriers( (1: 4) , : ) ;
Rx_carriers = Rx_carriers((5: 55), : ) ;
% --------------------------------------------- %
%                    信道估計                    %
% --------------------------------------------- %
Rx_training_symbols = Rx_training_symbols./ training_symbols;
Rx_training_symbols_deno = Rx_training_symbols.^2;
Rx_training_symbols_deno = Rx_training_symbols_deno(1,:)+Rx_training_symbols_deno(2,:)+Rx_training_symbols_deno(3,:)+Rx_training_symbols_deno(4,:) ;
Rx_training_symbols_nume = Rx_training_symbols(1, : ) +Rx_training_symbols(2, : ) + Rx_training_symbols(3, : ) +Rx_training_symbols(4, : ) ;
Rx_training_symbols_nume = conj(Rx_training_symbols_nume) ;
% % 取4個向量的導(dǎo)頻符號是為了進行平均優(yōu)化
% % 都是針對 “行向量”即單個的OFDM符號 進行操作
% % 原理：尋求1/H，對FFT之后的數(shù)據(jù)進行頻域補償
% % 1/H = conj(H)/H^2 because H^2 = H * conj(H) 
Rx_training_symbols = Rx_training_symbols_nume./Rx_training_symbols_deno;
Rx_training_symbols_2 = cat(1, Rx_training_symbols,Rx_training_symbols) ;
Rx_training_symbols_4 = cat(1, Rx_training_symbols_2,Rx_training_symbols_2) ;
Rx_training_symbols_8 = cat(1, Rx_training_symbols_4,Rx_training_symbols_4) ;
Rx_training_symbols_16 = cat(1, Rx_training_symbols_8, Rx_training_symbols_8) ;
Rx_training_symbols_32 = cat(1, Rx_training_symbols_16, Rx_training_symbols_16) ;
Rx_training_symbols_48 = cat(1, Rx_training_symbols_32, Rx_training_symbols_16) ;
Rx_training_symbols_50 = cat(1, Rx_training_symbols_48, Rx_training_symbols_2) ;
Rx_training_symbols = cat(1, Rx_training_symbols_50,Rx_training_symbols) ;
Rx_carriers = Rx_training_symbols.*Rx_carriers;
% 進行頻域單抽頭均衡
Rx_phase = angle(Rx_carriers)*(180/pi) ;
phase_negative = find(Rx_phase < 0) ;
%-------------------------------------------------------------------------
% 函數(shù)說明：
% %     FIND   Find indices of nonzero elements.
% %     I = FIND(X) returns the linear indices of the array X that are nonzero.
% %     X may be a logical expression. Use IND2SUB(I,SIZE(X)) to calculate
% %     multiple subscripts from the linear indices I.
%----------------------Test of Using "rem"---------------------------------
Rx_phase1 = Rx_phase; 
Rx_phase2 = Rx_phase;
Rx_phase1(phase_negative) = rem(Rx_phase1(phase_negative) + 360, 360) ;
Rx_phase2(phase_negative) = Rx_phase2(phase_negative) + 360 ;
if Rx_phase2(phase_negative) == Rx_phase1(phase_negative)
    fprintf('\n There is no need using rem in negative phase transition.\n')
else
    fprintf('\n We need to use rem in negative phase transition.\n')    
end
%-------------------------------------------------------------------------
Rx_phase(phase_negative) = rem(Rx_phase(phase_negative) + 360, 360) ;
% % 把負的相位轉(zhuǎn)化為正的相位
Rx_decoded_phase = diff(Rx_phase) ;
% %  這也是為什么要在前面加上初始相位的原因 
% % “Here a row vector of zeros is between training symbols and data symbols!!!”
%-------------------------------------------------------------------------
% 函數(shù)說明：
% %     DIFF Difference and approximate derivative.
% %     DIFF(X), for a vector X, is [X(2)-X(1)  X(3)-X(2) ... X(n)-X(n-1)].
% %     DIFF(X), for a matrix X, is the matrix of row differences,
% %        [X(2:n,:) - X(1:n-1,:)].
%------------------------Test Codes --------------------------------------
% %  a = [1 2 3; 4 5 6; 7 8 9; 10 11 12];
% %  b = a;
% %      for i = 2:4
% %        b(i,:) = b(i-1,:) + b(i,:);
% %      end
% %  c = diff(b);
%-----------------------Test Result --------------------------------------
% %  a =                  Modulating signal
% %      1     2     3
% %      4     5     6
% %      7     8     9
% %     10    11    12
% %  b =                  Modulated signal
% %      1     2     3
% %      5     7     9
% %     12    15    18
% %     22    26    30
% %  c =                  Recovered signal
% %      4     5     6
% %      7     8     9
% %     10    11    12
% % ----------------------------------------------------------------------------
% %   Name                   Size                    Bytes  Class
% %   Rx_phase              51x200                   81600  double array
% %   Rx_decoded_phase      50x200                   80000  double array
% % ---------------------------------------------------------------------------- 
phase_negative = find(Rx_decoded_phase < 0) ;
Rx_decoded_phase(phase_negative)= rem(Rx_decoded_phase(phase_negative) + 360, 360) ;
% % 再次把負的相位轉(zhuǎn)化為正的相位
% --------------------------------------------- %
%                   QDPSK解調(diào)                   %
% --------------------------------------------- %
base_phase = 360 /2^bits_per_symbol;
delta_phase = base_phase /2;
Rx_decoded_symbols = zeros(size(Rx_decoded_phase,1),size(Rx_decoded_phase,2)) ;
%
for i = 1: (2^bits_per_symbol - 1)
center_phase = base_phase*i;
plus_delta = center_phase + delta_phase;  % Decision threshold 1
minus_delta = center_phase - delta_phase; % Decision threshold 2
decoded = find((Rx_decoded_phase <= plus_delta)&(Rx_decoded_phase > minus_delta)) ;
Rx_decoded_symbols(decoded) = i;
end
% %  僅僅對三個區(qū)域進行判決
% %  剩下的區(qū)域就是零相位的空間了
% %  這個區(qū)域在定義解調(diào)矩陣時已經(jīng)定義為零
Rx_serial_symbols = reshape(Rx_decoded_symbols',1,size(Rx_decoded_symbols, 1)*size(Rx_decoded_symbols,2)) ;
for i = bits_per_symbol: -1: 1
if i ~= 1
Rx_binary_matrix(i, : ) = rem(Rx_serial_symbols, 2) ;
Rx_serial_symbols = floor(Rx_serial_symbols/2) ;
else
Rx_binary_matrix( i, : ) = Rx_serial_symbols;
end
end
% % Integer to binary
baseband_in = reshape(Rx_binary_matrix, 1,size(Rx_binary_matrix, 1)*size(Rx_binary_matrix, 2) ) ;
% --------------------------------------------- %
%                   誤碼率計算                   %
% --------------------------------------------- %
bit_errors = find(baseband_in ~= baseband_out) ;
% % find的結(jié)果 其每個元素為滿足邏輯條件的輸入向量的標號，其向量長度也就是收發(fā)不一樣的bit的個數(shù)
bit_error_count = size(bit_errors, 2) ;
total_bits = size( baseband_out, 2) ;
bit_error_rate = bit_error_count/ total_bits;
fprintf ( '%f \n',bit_error_rate) ;
% --------------------------------------------- %
%                   The END                     %
% --------------------------------------------- %