qam16_gain_pd.m

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

% qam16_Gain_PD.m
%
% Simulate 16QAM system with Gain PD
%
% 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: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.5;                              % 功率回退系数   
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);

%***************************************************************
%
%       Ideal Predistortion
%
%***************************************************************

% for n=1:length(PD_in_Env)
%     if (PD_in_Env(1,n)==0)
%         PD_in_Env(1,n)=0.000000001;
%     end;
% end;
%     PD_out_Env = (1-sqrt(1-PD_in_Env.^2))./PD_in_Env;   
%     
%     PD_out_Phase = (-pi/3)*(PD_out_Env.^2)./(1+PD_out_Env.^2) + PD_in_Phase;
% 
%     PD_out = PD_out_Env.*exp(j*PD_out_Phase);


%***************************************************************
%
%       增益预失真的查找表 ( LUT )
%
%***************************************************************

delta = 1/2^6;                                         % 查找表精度
t = [delta:delta:1];                                   % 表格长度
Fr = (1-sqrt(1-(t/2).^2))./((t/2).^2);                 % 存储的幅值
Fp = (-pi/3)*((Fr.*(t/2)).^2)./(1+(Fr.*(t/2)).^2);     % 存储的相位
F = Fr.*exp(j*Fp);                                     % 增益预失真的函数

   
%******************************************************************
%
%          增益预失真 ( Gain Based Predistortion )
%
%******************************************************************
 
num = fix(PD_in_Env/delta) + 1;                % 将信号对应到相应的表格
for n=1:length(num)
    FF(n) = F(num(n));
end;
PD_out = PD_in.*FF;


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

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(j*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;
     
    
  
   
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,'*-');

[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 ***********************************