ofdm_code.m

OFDM 编写源代码 matlab ofdm coder

% --------------------------------------------- %
% 摘录自《基于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
%-------------------------------------------------------------------------