AWGN信道条件下，基于16QAM+OFDM的误码率计算，并与理论值对比

1. 16QAM调制解调；
2. OFDM调制解调；
3. 不采用/采用卷积码；
4. AWGN信道

% an example of OFDM，Rayleigh+AWGN，理论值和仿真值的情况仿真
% Programmed by kkk, 2021/3/21
% OFDM参数
clear;
close all;
effc_symbols_num=1;                 % 用于存放有效数据的OFDM符号个数
carrier_num=64;                     % number of effective subcarriers in each OFDM Symbol
IFFT_bin_length=carrier_num;        % Length of FFT, which also means the subcarrier number in each OFDM Symbol
M_order=4;

PrefixRatio=1/4;                    % ratio of cyclic prefix, between 1/6~1/4
GI=PrefixRatio*carrier_num;         % The length of the guard interval
beta=1/32;                          % rcosdesign
% GIP=beta*(carrier_num+GI);        % The length of the Loop suffix
GIP=0;
EbN0=0:1:10;                        % 接收信噪比
SimTimes=1e2;
% 信道参数
PowerdB=[0 -8 -17 -21 -25];     % Channel tap power profile 'dB'
Delay=[0 3 5 6 8];              % Channel delay 'sample'
Power=10.^(PowerdB/10);         % Channel tap power profile 'linear scale'
Ntap=length(PowerdB);           % Chanel tap number
Lch=Delay(end)+1;               % 信道长度,是指一个信号传播时持续的时间，


% ChCodeType=0对应不采用信道编码；ChCodeType=1对应采用卷积编码
for ChCodeType=1:2
   for ee=1:length(EbN0)
       bers=0;
       for ss=1:SimTimes
           %************** Generate the baseband bits **********************
           if ChCodeType==1
               Xt=randi([0 1],effc_symbols_num*carrier_num*M_order,1);   % bit: integer vector
               encodedData=Xt;
               sys_symbols_num=effc_symbols_num;
           elseif ChCodeType==2                                          % 卷积编码(信道编码)
               Xt=randi([0 1],effc_symbols_num*carrier_num*M_order,1);   % bit: integer vector
               hConEnc = comm.ConvolutionalEncoder;
               encodedData = step(hConEnc, Xt);
               sys_symbols_num=length(encodedData)/(carrier_num*4);
           end
           Xt_after_Modulation=qam16(encodedData);
           Xt_after_Modulation_matrix=reshape(Xt_after_Modulation,carrier_num,sys_symbols_num);
           Xt_after_IFFT=ifft(Xt_after_Modulation_matrix);

           %************** Add the loop prefix/suffix **********************
           XX=zeros(IFFT_bin_length+GI+GIP,sys_symbols_num);
           for k=1:sys_symbols_num
               for i=1:IFFT_bin_length      % add the data
                   XX(i+GI,k)=Xt_after_IFFT(i,k); 
               end
               for i=1:GI                   % add the cyclic prefix
                   XX(i,k)=Xt_after_IFFT(i+IFFT_bin_length-GI,k);  
               end
               for i=1:GIP                  % add the cyclic suffix
                   XX(IFFT_bin_length+GI+i,k)=Xt_after_IFFT(i,k);  
               end    
           end           
           for i=1:sys_symbols_num
               windowed_Tx_data((IFFT_bin_length+GI)*(i-1)+1:(IFFT_bin_length+GI)*i)= XX(:,i);
           end
           
           %************** channel and AWGN **********************
           y= windowed_Tx_data;        
           Tx_signal_power=var(y);   %计算信号功率
           snr=EbN0(ee)+10*log10(M_order*(carrier_num/(carrier_num)));
           noise_mag = sqrt((10.^(-snr/10))*Tx_signal_power/2);
           Rx_data=y+noise_mag*(randn(size(y))+1i*randn(size(y)));

           %******************** OFDM Demoulation **********************
           Rx_data_matrix=zeros(IFFT_bin_length+GI+GIP,sys_symbols_num);           % 存放并行的接收数据
           for i=1:sys_symbols_num
               Rx_data_matrix(:,i)=Rx_data(1,(i-1)*(IFFT_bin_length+GI)+1:i*(IFFT_bin_length+GI)+GIP);
           end
           Rx_data_complex_matrix=Rx_data_matrix(GI+1:GI+IFFT_bin_length,:);
           Y1=fft(Rx_data_complex_matrix);
           Rx_phase=angle(Y1);                                                % 计算接收信号的相位特性
           Rx_mag=abs(Y1);                                                    % 计算接收信号的幅度特性
           [M,N] = pol2cart(Rx_phase,Rx_mag);                                          % 转换极坐标数据为直角坐标数据
           Rx_complex_carrier_matrix=complex(M,N);
           %******************** %16QAM 解调 **********************
           Rx_serial_complex_symbols=reshape(Rx_complex_carrier_matrix,size(Rx_complex_carrier_matrix,1)*size(Rx_complex_carrier_matrix,2),1);
           Xr=demoduqam16(Rx_serial_complex_symbols);          
           if ChCodeType==1
               errorRate = comm.ErrorRate;
               errors= step(errorRate, Xr, Xt);
               bers=bers+errors(1);
           elseif ChCodeType==2
               %******************** %卷积码解码 **********************
               hDec = comm.ViterbiDecoder('InputFormat','Hard');
               receivedBits = step(hDec, Xr);
               errorRate = comm.ErrorRate('ReceiveDelay',34);
               errors = step(errorRate, Xt, receivedBits);
               bers=bers+errors(1);       
           end      
       end
       ber(ChCodeType,ee)=bers/SimTimes;
   end
end

figure
% Analytic
EbN0dB=0:2:30;
ber_AWGN = ber_QAM(EbN0dB,2^M_order,'AWGN');
semilogy(EbN0dB,ber_AWGN,'r:');hold on;
% Simulation
semilogy(EbN0,ber(1,:),'bo-');hold on
semilogy(EbN0,ber(2,:),'ko-');hold on
grid on
legend('AWGN(Ana.)','AWGN (Sim.)', 'AWGN with Convolutional Code(Sim.)');
xlabel('EbN0(dB)');ylabel('SER');
axis([EbN0(1) EbN0(end) 1e-5 1]);

%========调用已有的子函数===============%
function [complex_qam_data]=qam16(bitdata)
% (input)  bitdata--bitstream
% (output) complex_qam_data--16QAM complex data
X1=reshape(bitdata,4,length(bitdata)/4)';
d=1;
for i=1:length(bitdata)/4
   for j=1:4
       X1(i,j)=X1(i,j)*(2^(4-j));
   end
   source(i,1)=1+sum(X1(i,:));
end
mapping=[-3*d   3*d;      -d	3*d;     d	3*d;      3*d	3*d;...
         -3*d	d;        -d	d;       d	d;        3*d	d;...
         -3*d	-d;       -d	-d;      d	-d;       3*d	-d;...
         -3*d	-3*d;     -d	-3*d;    d	-3*d;     3*d	-3*d];
for i=1:length(bitdata)/4
    qam_data(i,:)=mapping(source(i),:);
end
complex_qam_data=complex(qam_data(:,1),qam_data(:,2));

% 16QAM信号的解调子程序
function [demodu_bit_symble]=demoduqam16(Rx_serial_complex_symbols)
% 输入参数：Rx_serial_complex
% 输出参数：二进制数据码流
complex_symbols=reshape(Rx_serial_complex_symbols,length(Rx_serial_complex_symbols),1);
d=1;

mapping=[-3*d   3*d;      -d	3*d;     d	3*d;      3*d	3*d;...
         -3*d	d;        -d	d;       d	d;        3*d	d;...
         -3*d	-d;       -d	-d;      d	-d;       3*d	-d;...
         -3*d	-3*d;     -d	-3*d;    d	-3*d;     3*d	-3*d];
complex_mapping=complex(mapping(:,1),mapping(:,2));
for i=1:length(Rx_serial_complex_symbols)
   for j=1:16
      metrics(j)=abs(complex_symbols(i,1)-complex_mapping(j,1)); 
   end
   [min_metric decode_symble(i)]=min(metrics);
end
decode_bit_symble=de2bi((decode_symble-1)','left-msb');     % 十进制数转二进制数
demodu_bit_symble=reshape(decode_bit_symble',length(Rx_serial_complex_symbols)*4,1);


function ber=ber_QAM(EbN0dB,M,AWGN_or_Rayleigh)
% Find ananlytical BER of Mary QAM in AWGN or Rayleigh channel
% EbN0dB=EbN0dB: Energy per bit-to-noise power[dB] for AWGN channel
%       =rdB: Average SNR(2*sigma Eb/N0)[dB] for Rayleigh channel
% M = Modulation order (Constellation size)  

%MIMO-OFDM Wireless Communications with MATLAB㈢   Yong Soo Cho, Jaekwon Kim, Won Young Yang and Chung G. Kang
%?2010 John Wiley & Sons (Asia) Pte Ltd

N= length(EbN0dB);  sqM= sqrt(M); 
a= 2*(1-power(sqM,-1))/log2(sqM);  b= 6*log2(sqM)/(M-1);
if nargin<3,  AWGN_or_Rayleigh='AWGN';  end
if lower(AWGN_or_Rayleigh(1))=='a', ber = a*Q(sqrt(b*10.^(EbN0dB/10)));
 else   rn= b*10.^(EbN0dB/10)/2; ber = 0.5*a*(1-sqrt(rn./(rn+1)));
end

function y=Q(x)
% co-error function: 1/sqrt(2*pi) * int_x^inf exp(-t^2/2) dt.
y=erfc(x/sqrt(2))/2;

说明：AWGN(Ana.)为理论值；
AWGN(Sim.)为仿真值;
AWGN with Convolutional Code(Sim.)为采用卷积编码的仿真值