qam16_twta.m

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

% qam16_TWTA.m
%
% Simulate 16 QAM system with TWTA
%
% 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);
    
%***************************************************************
%
%       预失真之前的归一化和功率回退
%
%***************************************************************
IBO_dB =4;                              % 功率回退系数   
nf_ibo = 10^(-IBO_dB/10);                % 功率回退复系数

Ht_out = ich2 +i * qch2;
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);


%******************************************************************
%
%         行波管放大器 ( TWTA )
%
%*******************************************************************

PA_out_r = Ar*PD_in_Env./(1+Br*PD_in_Env.^2);                                  % 幅度非线性放大
PA_out_p = Ap*PD_in_Env.^2./(1+Bp*PD_in_Env.^2) + PD_in_Phase;               % 相位非线性放大

PA_out = PA_out_r.*exp(j*PA_out_p)/nf_ibo*nf/2;     % 去归一化和去功率回退
PA_out_i = real(PA_out);
PA_out_q = imag(PA_out);
     
for ebn0 = 1:length(SNR_dB)+1
     
%******************** START CALCULATION *****************************

nloop = 100;                       % 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;
     
     ber(ebn0) = noe/nod;
  
   
end  % for iii = 1:nloop

%******************** Output result *********************************

 
 %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-');

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

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


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

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