path3.m

用MATLAB对ofdm系统的性能仿真

% 1. 该程序进行了简单的LMS信道估计，没有加入与MMSE等其他信道估计算法的比较；
% 2. 仿真条件为系统处于理想同步情况下。
%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 All-zero CP  
% 调制方式            QDPSK
% 多径信道数          2、3、4（缺省）
% 信道最大时延        7 (单位数据符号)
% 仿真条件            收发之间严格同步 
%SNR = input('SNR =') ;
%fprintf('SNR= %d \n',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));   % 坐标： (1 to 200) + 156 ,  157 -- 356
conjugate_carriers = IFFT_bin_length - carriers + 2;    % 坐标 ：1024 - (157:356) + 2 = 1026 - (157:356) = (869:670)     
% 构造共轭时间－载波矩阵，以便应用所谓的RCC，Reduced Computational Complexity算法，即ifft之后结果为实数 
% Define the conjugate time-carrier matrix
% 也可以用flipdim函数构造对称共轭矩阵
% --------------------------------------------- %
%                   信号发射                     %
% --------------------------------------------- %
out = rand(1,baseband_out_length);
baseband_out1 = round(out) ;
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) ;
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'方式
%-------------------------------------------------------------------------
%  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');
%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) ;
% 差分调制 
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’.  
complex_carrier_matrix = complex(X, Y) ;
%－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－－
% 添加训练序列
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
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) ;
%-------------------------------------------------------------------------
%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';
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
ofdm_modulation = reshape(windowed_time_wave_matrix',1, IFFT_bin_length*(4 + symbols_per_carrier + 1) ) ;
ofdm_modulation_tmp = windowed_time_wave_matrix.';
ofdm_modulation_test = ofdm_modulation_tmp(:)';
%if (ofdm_modulation_test == ofdm_modulation)
%  fprintf('ofdm_modulation_test == ofdm_modulation \n\n\n');
%else
%  fprintf('ofdm_modulation_test ~= ofdm_modulation \n\n\n');
%end 
% We get the result "ofdm_modulation_test == ofdm_modulation" .
%-------------------------------------------------------------------------
Tx_data = ofdm_modulation;
% --------------------------------------------- %
%                   信道模拟                     %
% --------------------------------------------- %
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
Tx_data = Tx_data + copy1 + copy2 + copy3;
% 3 multi-paths
Tx_signal_power = var(Tx_data);
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;
Rx_Data = Tx_data + noise;
% --------------------------------------------- %
%                  信号接收                      %
% --------------------------------------------- %
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个向量的导频符号是为了进行平均优化
% % 都是针对 “行向量”即单个的OFDM符号 进行操作
% % 原理：寻求1/H，对FFT之后的数据进行频域补偿
% % 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) ;
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) ;
% % 把负的相位转化为正的相位
Rx_decoded_phase = diff(Rx_phase) ;
% %  这也是为什么要在前面加上初始相位的原因 
% % “Here a row vector of zeros is between training symbols and data symbols!!!”
phase_negative = find(Rx_decoded_phase < 0) ;
Rx_decoded_phase(phase_negative)= rem(Rx_decoded_phase(phase_negative) + 360, 360) ;
% % 再次把负的相位转化为正的相位
% --------------------------------------------- %
%                   QDPSK解调                   %
% --------------------------------------------- %
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
% %  仅仅对三个区域进行判决
% %  剩下的区域就是零相位的空间了
% %  这个区域在定义解调矩阵时已经定义为零
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的结果 其每个元素为满足逻辑条件的输入向量的标号，其向量长度也就是收发不一样的bit的个数
bit_error_count = size(bit_errors, 2) ;
total_bits = size( baseband_out, 2) ;
bit_error_rate = bit_error_count/ total_bits;
fprintf ( 'BER=%f (3 paths)\n',bit_error_rate) ;
% --------------------------------------------- %
%                   The END                     %
% --------------------------------------------- %