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