Matlab簡單的OFDM仿真,信道估計,有BER曲線

1、clear all;close all;fprintf( 'n OFDM仿真n n') ;% - %                   參數(shù)定義                     % - %IFFT_bin_length = 1024;carrier_count   = 200;bits_per

2、_symbol = 2;symbols_per_carrier = 50;% 子載波數(shù)            200% 位數(shù)/ 符號          2% 符號數(shù)/ 載波        50% 訓(xùn)練符號數(shù)          10% 循環(huán)前綴長度        T/4（作者注明）

3、  All-zero CP  % 調(diào)制方式            QDPSK% 多徑信道數(shù)          2、3、4（缺省）% 信道最大時延        7 (單位數(shù)據(jù)符號)% 仿真條件            收發(fā)之間嚴格同步 %SNR=input('SNR=

4、9;);    % 輸入信噪比參數(shù)SNR=3:14;%定義信噪比范圍BER=zeros(1,length(SNR);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);   % 坐標(biāo)： (1 to 200) + 156 ,  157 - 356c

5、onjugate_carriers=IFFT_bin_length-carriers+2;  % 坐標(biāo) ：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)造對稱共軛矩陣% - %         

6、0;         信號發(fā)射                     % - %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);% 四種生成發(fā)送的二進制序列的方法，任取一種產(chǎn)生要發(fā)送的二進制序列%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% 驗證四種生成發(fā)送的二進

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個比特轉(zhuǎn)化為整數(shù) 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');% 函數(shù)說明：% 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調(diào)制                   % - %carrier_matrix = zeros(1,carrier_count); carrier_matrix;  % 添加一個差分調(diào)制的初始相位，為0for i = 2symbols_per_carrier + 1)    carrier_matrix(i, = rem(carrier_m

15、atrix(i, + carrier_matrix (i-1, 2bits_per_symbol) ;  % 差分調(diào)制 endcarrier_matrix = carrier_matrix*(2*pi)/(2bits_per_symbol) ;  % 產(chǎn)生差分相位X, Y=pol2cart(carrier_matrix, ones(size(carrier_matrix,1),size(carrier_matrix,2); % 由極坐標(biāo)向復(fù)數(shù)坐標(biāo)轉(zhuǎn)化 第一參數(shù)為相位 第二參數(shù)為幅度% Carrier_matrix contains all the pha

16、se information and all the amplitudes are the same1. complex_carrier_matrix = complex(X, Y) ;% 添加訓(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 -

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) ; % 訓(xùn)練序列與數(shù)據(jù)合并 % 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結(jié)果都是實數(shù)% 由此可以看出，只是取了IFFT之后載波上的點，并未進行CP的復(fù)制和添加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個向量的導(dǎo)頻符號是為了進行平均優(yōu)化% 都是針對 “行向量”即單個的OFDM符號 進行操作% 原理：

33、尋求1/H，對FFT之后的數(shù)據(jù)進行頻域補償% 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; % 把負的相位轉(zhuǎn)化為正的相位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、;  % 再次把負的相位轉(zhuǎn)化為正的相位% - %                   QDPSK解調(diào)                   % - % 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%  僅僅對三個區(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 =