qam16_secant.m

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

% qam16_secant.m
%
% Simulation 16QAM system with secant 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 = 2^10;                     % Number of symbols htat simulates in
                               % each loop
IPOINT = 8;                    % Number of oversamples

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

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 )
%
%******************************************************************
 
delta = 1/2^6;
F = zeros(1,1/delta);    % 将信号对应到相应的表格

%******************************************************************
%
%          自适应 ( secant method )
%
%******************************************************************

tic;
for n = 1:length(PD_in)
aa = n
num(n) = fix(PD_in_Env(n)/delta) + 1;

 F0 = 0;
 F1 = 1;
  
 while abs(F1-F0) > 0.1
    
      PD0(n) = PD_in(n) * F0;
      PD1(n) = PD_in(n) * F1;                   

      PA0a(n) = Ar*abs(PD0(n))/(1+Br*abs(PD0(n))^2);                                  % 幅度非线性放大
      PA0p(n) = Ap*abs(PD0(n))^2/(1+Bp*abs(PD0(n))^2) + angle(PD0(n));               % 相位非线性放大
      PA0(n) = PA0a(n)*exp(i*PA0p(n))/nf_ibo*nf;    % 去归一化和去功率回退

      PA1a(n) = Ar*abs(PD1(n))/(1+Br*abs(PD1(n))^2);                                  % 幅度非线性放大
      PA1p(n) = Ap*abs(PD1(n))^2/(1+Bp*abs(PD1(n))^2) + angle(PD1(n));               % 相位非线性放大
      PA1(n) = PA1a(n)*exp(i*PA1p(n))/nf_ibo*nf;     % 去归一化和去功率回退
          
      e0 =  PA0(n) - Ht_out(n);
      e1 =  PA1(n) - Ht_out(n);
  
      if (e1-e0) == 0
            fenmu = 0.1;
      else 
             fenmu = e1 - e0;   
      end
  
      F2 = (F0*e1 - F1*e0)/fenmu;
      F0 = F1;
      F1 = F2;

end  % while

FF(n) = F2;

end;
toc;

 PA_out_i = real(PA1);
 PA_out_q = imag(PA1);

     
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,'v-');

[pxx1,f1] = pwelch(Ht_out,512);                 
pxxdB1 = 10 * log10(pxx1 / max(pxx1));

[pxx2,f2] = pwelch(PA1,512);                 
pxxdB2 = 10 * log10(pxx2 / max(pxx2));
figure(4);                                          % OFDM调制、解调信号功率谱密度
plot(f2,pxxdB2,'y');

eyediagram((ich+qch*i),2);
eyediagram((ich5(1:500)+qch5(1:500)*i),2)

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