%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%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