qam16_discretenewton.m

上传的是一个16Qam预失真处理的程序

% qam16_DiscreteNewton.m
%
% Simulate16QAM with Discrete Newton's method adaptive algorithm
%
% Programmed by linxiaochen  
%

%******************** Preparatin part *******************************

sr = 256000.0;                 % Symbol rate
m1 = 4;                        % m1:Number of modulation levels
                               % (BPSK:m1=1, QPSK:m1=2, 16QAM:m1=4)
br = sr .* m1;                 % Bit rate
nd = 1024;                     % Number of symbols htat simulates in
                               % each loop
IPOINT = 8;                    % Number of oversamples

SNR_dB = 1:12;                           % 仿真信噪比范围
SNR1_dB = 0:0.1:12;  

Ar = 2.0;                                % TWTA的Dong-Seog Han参数
Br = 1;
Ap = pi/3;
Bp = 1;

%******************** Filter initialization *************************

irfn = 21;                       % Number of taps
alfs = 0.5;                      % Rolloff factor
[xh] = hrollfcoef(irfn,IPOINT,sr,alfs,1);
    % Transmitter filter coefficients
[xh2] = hrollfcoef(irfn,IPOINT,sr,alfs,0);
    % Receiver filter coefficients
 
%******************** Data generation *******************************
     
      data1 = rand(1,nd*m1) > 0.5;    % rand:built in function

%******************** 16QAM Modulation ******************************

     [ich,qch] = qammod(data1,1,nd,m1);
     figure(1);
      plot(ich,qch,'*');
     [ich1,qch1] = compoversamp(ich,qch,length(ich),IPOINT);
     [ich2,qch2] = compconv(ich1,qch1,xh);
 
     Ht_out = ich2 +i * qch2;
     
%***************************************************************
%
%       预失真之前的归一化和功率回退
%
%***************************************************************
IBO_dB = 4.5;                              % 功率回退系数   
nf_ibo = 10^(-IBO_dB/10);                % 功率回退复系数

nf = sqrt(0.5*mean(abs(Ht_out).^2));                    % 归一化系数

PD_in = nf_ibo*Ht_out/nf;                             % 归一化和功率回退

PD_in_Env = abs(PD_in);
PD_in_Phase = angle(PD_in);

%******************************************************************
%
%          增益预失真 ( Gain Based Predistortion )
%
%******************************************************************
 
               
k = zeros(1,length(PD_in_Env));
tic;
for n = 1:length(Ht_out)
aa = n

 F0 = [0;0];
 F1 = [1;1];


 while sum(abs(F1-F0)) > 0.1

    PD_out_01 = PD_in(n)*F0(1).*exp(i*F1(2));
    PD_out_10 = PD_in(n)*F1(1).*exp(i*F0(2));
    PD_out_11 = PD_in(n)*F1(1).*exp(i*F1(2));

    
    PA_out_01_r = Ar*abs(PD_out_01)./(1+Br*abs(PD_out_01).^2);                                  % 幅度非线性放大
    PA_out_01_p = Ap*abs(PD_out_01).^2./(1+Bp*abs(PD_out_01).^2) + angle(PD_out_01);               % 相位非线性放大 
    
    PA_out_10_r = Ar*abs(PD_out_10)./(1+Br*abs(PD_out_10).^2);                                  % 幅度非线性放大
    PA_out_10_p = Ap*abs(PD_out_10).^2./(1+Bp*abs(PD_out_10).^2) + angle(PD_out_10);               % 相位非线性放大 
     
    PA_out_11_r = Ar*abs(PD_out_11)./(1+Br*abs(PD_out_11).^2);                                  % 幅度非线性放大
    PA_out_11_p = Ap*abs(PD_out_11).^2./(1+Bp*abs(PD_out_11).^2) + angle(PD_out_11);               % 相位非线性放大
            
   M00 = (PA_out_11_r - PA_out_01_r)/(F1(1) - F0(1));
   M01 = (PA_out_11_r - PA_out_10_r)/(F1(2) - F0(2));
   M10 = (PA_out_11_p - PA_out_01_p)/(F1(1) - F0(1));
   M11 = (PA_out_11_p - PA_out_10_p)/(F1(2) - F0(2));
            

    e = [PA_out_11_r-abs(PD_in(n));PA_out_11_p-angle(PD_in(n))];
    
     M_ni = [M11,-M01;-M10,M00]/(M00*M11 - M01*M10);

     F2 = F1 - M_ni*e;
     F0 = F1;
     F1 = F2;
 

end  % while

PD_out(n) = PD_in(n)*F0(1)*exp(i*F0(2));

end;
toc;

PA_out_r = Ar*abs(PD_out)./(1+Br*abs(PD_out).^2);                                  % 幅度非线性放大
PA_out_p = Ap*abs(PD_out).^2./(1+Bp*abs(PD_out).^2) + angle(PD_out);               % 相位非线性放大
         
PA_out =PA_out_r.*exp(i*PA_out_p)/nf_ibo*nf;
PA_out_i = real(PA_out);
PA_out_q = imag(PA_out);

     
for ebn0 = 1:length(SNR_dB)+1
     
%******************** START CALCULATION *****************************

nloop = 10;                       % Number of simulation loops
noe = 0;                           % Number of error data
nod = 0;                           % Number of transmitted data

  for iii = 1:nloop     
%******************** Attenuation Calculation ***********************

      spow = sum(PA_out_i.*PA_out_i+PA_out_q.*PA_out_q)/nd;
          % sum:built in function
      attn = 0.5*spow*sr/br*10.^(-(ebn0-1)/10);
      attn = sqrt(attn);
          % sqrt:built in function
     
%************** Add White Gaussian Noise (AWGN) *********************

     [ich3,qch3] = comb(PA_out_i,PA_out_q,attn);
          % add white gaussian noise
     [ich4,qch4] = compconv(ich3,qch3,xh2);
     
     sampl = irfn*IPOINT+1;
     ich5 = ich4(sampl:IPOINT:length(ich4));
     qch5 = qch4(sampl:IPOINT:length(qch4));
     
     ich6 = ich5(1:1000);
     qch6 = qch5(1:1000);
     figure(2);
     plot(ich6,qch6,'*');
     
%******************** 16QAM Demodulation ****************************
     
     [demodata] = qamdemod(ich5,qch5,1,nd,m1);
      
%******************** Bit Error Rate (BER) **************************
  
     noe2 = sum(abs(data1-demodata));
     nod2 = length(data1);
     noe = noe + noe2;
     nod = nod + nod2;
     
     
    % fprintf:built in function
   
end  % for iii = 1:nloop

%******************** Output result *********************************
 ber(ebn0) = noe/nod;
 
end

t1 = [0:0.1:12];
tt1=exp(t1*log(10)/10);  
B1 = 3/8.*erfc(sqrt(2/5.*tt1))-9/64.*erfc(sqrt(2/5.*tt1)).*erfc(sqrt(2/5.*tt1));
t11=[0:length(ber)-1];
figure(3);
semilogy(t1,B1,t11,ber,'o-');

%******************** end of file ***********************************