%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%main 主程序
clear;
close all;
numsymb=256*2;
M=16; % mapping point
N=10; % sampling frequency
fd=512;
fs=N*fd;
msg_orig=randsrc(numsymb,1,[0:M-1],4321);
msg_tx=modmap(msg_orig,fd,fd,'qam',M); %constellation mapping
x=complex(msg_tx(:,1),msg_tx(:,2));
%raised cosine filter
R=0.5;delay=6;
yf=rcosine(fd,fs,'fir',R,delay);
x_f=rcosflt(x,fd,fs,'filter',yf);
msg_x=x_f(1+delay*N:end-delay*N);
% figure(1)
% plot(1:10:200,real(x(1:20)),'o')
% hold on;
% plot(1:1:200,real(msg_x(1:200)),'r-')
% msg_x=x_f;
msg_x=msg_x/max(abs(msg_x))*1;
ww=PA_model(msg_x);
zz=PD_deal(msg_x);
yy=PA_model(zz);
SNR=60;
yy=awgn(yy,SNR,'measured',1234,'dB');
% h1=[];
% h1=scatterplot(yy,N,0,'r.',h1);
[px1,ff]=pwelch(yy,1024*2,512*2,1024*2,fs);
pxx1=fftshift(px1);
[px2,ff]=pwelch(ww,1024*2,512*2,1024*2,fs);
pxx2=fftshift(px2);
[px3,ff]=pwelch(msg_x,1024*2,512*2,1024*2,fs);
pxx3=fftshift(px3);
ff=(0:length(px1)-1)*fs/fd/length(px1)-fs/fd/2;
pxx1_db=10*log10(abs(pxx1));
pxx2_db=10*log10(abs(pxx2));
pxx3_db=10*log10(abs(pxx3));
p1=decimate(pxx1_db,1);
p2=decimate(pxx2_db,1);
p3=decimate(pxx3_db,1);
f=decimate(ff,1);
figure(4);
plot(f,p1,'r');grid,hold on;
plot(f,p2,'b');grid,hold on;
%plot(f,p3,'g');grid,hold on;
xlabel('Normalized frequency (f-fc)/f1');
ylabel('Output power spectrum (dB)');
[ACPR1, AltPR1]=acpr_comp(ff,p1)%预失真
[ACPR2, AltPR2]=acpr_comp(ff,p2)%PA输入输出
[ACPR3, AltPR3]=acpr_comp(ff,p3)%
%%%%%%%%%%%%%%%%% PA_model.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%PA模型
function [pa_out]=PA_model(msg)
%msg=0:0.1:1;
r=abs(msg);
a1=2.1587;a2=1.1517;
mag_out=a1*r./(1+a2.*r.*r);
%for n=1:length(msg)
% if mag_out(n)>1
% mag_out(n)=1;
% end
%end
b1=4.033;b2=9.1040;
delta_phase=b1.*r.*r./(1+b2.*r.*r);
phase=delta_phase+angle(msg);
pa_out=mag_out.*exp(i*phase);
%%%%%%%%%%%%%%%%%%% PD_model.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%PD预失真模型
function [pd_para1,pd_para2]=(msg)
a=0.05;
v=[1 1 1];
t=[1 1 1];
r1=abs(msg);
num=length(msg);
r0=r1.^0;
r3=r1.^3;
r5=r1.^5;
r=[r1.';r3.';r5.'];
s=[r0.';r1.';r3.'];
for i=1:num
%abs;
m=v*r(:,i);
n=2.1587.*m./(1+1.1517.*m.*m);
p=r1(i,:)-n;
q=a.*p.*r(:,i);
v=v+q.';
%angle;
d=4.033.*m.*m./(1+9.1040.*m.*m);
f=t*s(:,i);
e=-d-f;
g=a.*e.*s(:,i);
t=t+g.';
end
pd_para1=v;
pd_para2=t;
%%%%%%%%%%%% PD_deal.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%PD预失真处理
function [pd_out]=PD_deal(msg)
r1=abs(msg);
r0=r1.^0;
r3=r1.^3;
r5=r1.^5;
r=[r1.';r3.';r5.'].';
s=[r0.';r1.';r3.'].';
[p1, p2]=PD_model(msg);
%abs;
m=r*p1.';
%angle;
f=s*p2.';
phase=f+angle(msg);
pd_out=m.*exp(i*phase);
%%%%%%%%%%%% acpr_comp.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%acpr计算
function [ACPR, AltPR]=acpr_comp(ff,p
