简单的OFDM仿真信道估计有BER曲线

1、clear all;close all;fprintf( 'n OFDM仿真n n') ;% - %                   参数定义                     % - %IFFT_bin_length = 1024;carrier_count   = 200;bits_per

2、_symbol = 2;symbols_per_carrier = 50;% 子载波数            200% 位数/ 符号          2% 符号数/ 载波        50% 训练符号数          10% 循环前缀长度        T/4（作者注明）

3、  All-zero CP  % 调制方式            QDPSK% 多径信道数          2、3、4（缺省）% 信道最大时延        7 (单位数据符号)% 仿真条件            收发之间严格同步 %SNR=input('SNR=

4、9;);    % 输入信噪比参数SNR=3:14;%定义信噪比范围BER=zeros(1,length(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 - 356c

5、onjugate_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函数构造对称共轭矩阵% - %         

6、0;         信号发射                     % - %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_ou

7、t_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% 验证四种生成发送的二进

8、制序列的方法baseband_out=round( rand(1,baseband_out_length);convert_matrix = reshape(baseband_out,bits_per_symbol,length(baseband_out)/bits_per_symbol);for k = 1length(baseband_out)/bits_per_symbol),  modulo_baseband(k) = 0; for i = 1:bits_per_symbol     modulo_baseband(k) =

9、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)

10、/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

11、= 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 thebinary% input is Right-MSB; the first element in B represents the least significant bit.%if (mod

12、ulo_baseband = Test_convert_matrix1') %   fprintf('modulo_baseband = Test_convert_matrix1 nnn');%else if (modulo_baseband = Test_convert_matrix2')     %    fprintf('modulo_baseband = Test_convert_matrix2 nnn');%    else% 

13、   fprintf('modulo_baseband = any Test_convert_matrix nnn'); %    end%end% we get the result "modulo_baseband = Test_convert_matrix1".%-carrier_matrix = reshape(modulo_baseband,carrier_count,symbols_per_carrier)'% 生成时间载波矩阵% - %       

14、60;          QDPSK调制                   % - %carrier_matrix = zeros(1,carrier_count); carrier_matrix;  % 添加一个差分调制的初始相位，为0for i = 2symbols_per_carrier + 1)    carrier_matrix(i, = rem(carrier_m

15、atrix(i, + carrier_matrix (i-1, 2bits_per_symbol) ;  % 差分调制 endcarrier_matrix = carrier_matrix*(2*pi)/(2bits_per_symbol) ;  % 产生差分相位X, Y=pol2cart(carrier_matrix, ones(size(carrier_matrix,1),size(carrier_matrix,2); % 由极坐标向复数坐标转化 第一参数为相位 第二参数为幅度% Carrier_matrix contains all the pha

16、se information and all the amplitudes are the same1. 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 -

17、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 -

18、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 rowtraining_symbols = cat(1, training_symbols, training_symbols) ;  training_symbols = cat(1, training_symbols, training_symbols) ; % Production of 4 rows of training_symb

19、olscomplex_carrier_matrix = cat(1, training_symbols, complex_carrier_matrix) ; % 训练序列与数据合并 % block-type pilot symbolsIFFT_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

20、% 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_car

21、rier,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_carrie

22、r%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&#

23、160;  windowed_time_wave_matrix( i, : ) = real(time_wave_matrix( i, : ) ;end% get the real part of the result of IFFT% 这一步可以省略，因为IFFT结果都是实数% 由此可以看出，只是取了IFFT之后载波上的点，并未进行CP的复制和添加endofdm_modulation = reshape(windowed_time_wave_matrix',1, IFFT_bin_length*(4 + symbols_per_carrier + 1) ) ;% P2S o

24、peration%-%   Test by lavabin%   Another way of matrix transition%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 nnn');%els

25、e%fprintf('ofdm_modulation_test = ofdm_modulation nnn');%end % We get the result "ofdm_modulation_test = ofdm_modulation" .%-Tx_data=ofdm_modulation;% - %                   信道模拟          

26、           % - %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) ;endcopy2 = zeros(size(Tx_data) ) ;for i = 1 +

27、d2: length( Tx_data)copy2(i) = a2*Tx_data( i - d2) ;endcopy3 = zeros(size(Tx_data) ) ;for i = 1 + d3: length(Tx_data)copy3(i) = a3*Tx_data ( i - d3) ;endcopy4 = zeros(size(Tx_data) ) ;for i = 1 + d4: length( Tx_data)copy4(i) = a4*Tx_data(i - d4) ;endTx_data = Tx_data + copy1 + copy2 + copy3 + copy4;

28、 % 4 multi-pathsTx_signal_power = var(Tx_data);for idx=1:length(SNR)%monte carlo 仿真模拟    linear_SNR = 10( SNR(idx) /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;% - %&

29、#160;                 信号接收                      % - % Rx_Data_matrix = reshape(Rx_Data, IFFT_bin_length, 4 + symbols_per_carrier + 1) ;Rx_spectrum = fft(Rx_Data_matrix) ; %

30、60; Suppose precise synchronazition between Tx and RxRx_carriers = Rx_spectrum( carriers, : )'Rx_training_symbols = Rx_carriers( (1: 4) , : ) ;Rx_carriers = Rx_carriers(5: 55), : ) ;% - %                    信道估计 

31、60;                 % - %   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_traini

32、ng_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符号 进行操作% 原理：

33、寻求1/H，对FFT之后的数据进行频域补偿% 1/H = conj(H)/H2 because H2 = H * conj(H) Rx_training_symbols = Rx_training_symbols_nume./Rx_training_symbols_deno;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_symb

34、ols_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

35、_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_phas

36、e = angle(Rx_carriers)*(180/pi) ;phase_negative = find(Rx_phase < 0) ;%-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(ph

37、ase_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) ;&

38、#160; % 把负的相位转化为正的相位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)

39、;  % 再次把负的相位转化为正的相位% - %                   QDPSK解调                   % - % base_phase = 360 /2bits_per_symbol;delta_phase = base_phase /2;Rx_decoded_symbols = zeros(size(

40、Rx_decoded_phase,1),size(Rx_decoded_phase,2) ;for i = 1: (2bits_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

41、 = 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 =