clear;
clc;
carr=1024; %子载波个数
sym_c=10; %每个子载波含有的OFDM符号数
bit_sym=2; %每个符号含有的比特数，QPSK调制
IFFT_n=256; %IFFT点数
r=1/10; %保护间隔和OFDM数据的比例；
%产生信号
sum=carr*sym_c*bit_sym;
colume=sum/(2*carr);
signal=rand(1,sum)>0.5; %初始信号
%QPSK调制,QPSK_sig里面存放的是调制后的信号，数目sumQ
sumQ=sum/2;%OFDM符号总数
imag=sqrt(-1); % 虚部 j
QPSK=[-1+imag,-1-imag,1+imag,1-imag]; %创建QPSK 映射表
SIGNAL=zeros(1,sumQ); %计算并存放调制前的十进制数据(将二进制转为十进制数据，第20行)
QPSK_sig=zeros(1,sumQ); %存放调制后的QPSK信号
for n=1:sumQ
SIGNAL(n)=signal(2*n-1)*2+signal(2*n); %将二进制换算成十进制
end
for i=1:sumQ
if SIGNAL(i)==0;
QPSK_sig(i)=QPSK(1);
elseif SIGNAL(i)==1;
QPSK_sig(i)=QPSK(2);
elseif SIGNAL(i)==2;
QPSK_sig(i)=QPSK(3);
elseif SIGNAL(i)==3;
QPSK_sig(i)=QPSK(4);
end
end
%串/并转换 计算第i个载波上面的信号to_par（i，：）
colume=sumQ/carr;%每个子载波的符号数
for i=1:carr % carr载波个数
for j=1:colume;
to_par(i,j)=QPSK_sig(carr*(j-1)+i);
end
end
colume=sumQ/carr;
% % to_par=reshape(QPSK_sig,carr,colume);
%每个子载波上进行 IFFT变换 （调制后的QPSK信号进行IFFT）
for j=1:colume
y(:,j)=ifft(to_par(:,j)); %每一列进行一次IFFT
end
% % y=ifft(to_par);
yr=real(y); %实部
yi=(y-yr)*sqrt(-1)*(-1);%虚部
y_=yr+yi*sqrt(-1);
y_series=reshape(y_,1,(colume)*carr);
y1=y_(:,1);
y2=abs(y1);%采样点幅值
y2_power=y2.^2;%采样点功率
y2_power_mean=mean(y2_power);
sigma=sqrt(y2_power_mean);
y_PAPR=10*log10(y2_power./sigma^2);
PAPR_max=max(y_PAPR);
PAPR_mean=mean(y_PAPR);
%*******************************************进行TR消峰处理***********************
L=32;%预留子载波数目
u=0.3;%设置步进参数
for i=1:carr; %生成IFFT变化矩阵
for j=1:carr;
Q(i,j)=exp(sqrt(-1)*2*pi*(i-1)*(j-1)/carr)/sqrt(carr);
end
end
Ql=Q(:,carr-L+1:carr); %矩阵Q尖
x(:,1)=y1;%步骤一：算法初始化
Imax=16; %迭代次数
i=0;
A=y2_power_mean;%门限值
U=zeros(carr,1);
flag=1;
while((i<Imax)&&(flag==1))
flag=0;
for n=1:carr
if abs(x(n))>A
q_nrow=Ql(n,:);
U=U+(x(n)-A*exp(sqrt(-1)*phase(x(n))))*Ql*q_nrow';
flag=1;
else
U=U+zeros(carr,1);
end
end
x=x-u*U;
i=i+1;
end
x_power=abs(x).^2;
x_power_mean=mean(x_power);
x_sigma=sqrt(x_power_mean);
x_PAPR=10*log10(x_power./x_sigma^2);
subplot(2,1,1);
plot(y2_power,'b-');
xlim([0,1024]);
ylim([0,0.02]);
grid on;
subplot(2,1,2);
plot(x_power,'r-');
xlim([0,1024]);
ylim([0,0.02]);
grid on;
figure;
%************************************统计CCDF特性*************************************
PAPR0=2:0.5:9;
for k=1:length(PAPR0)
count(k)=0;
x_count(k)=0;
end
for l=1:carr;
for k=1:length(PAPR0)
if y_PAPR(l)>PAPR0(k)
count(k)=count(k)+1;
end
if x_PAPR(l)>PAPR0(k)
x_count(k)=x_count(k)+1;
end
end
end
CCDF=count/carr;
x_CCDF=x_count/carr;
subplot;
semilogy(PAPR0,CCDF,'.-k');
hold on;
semilogy(PAPR0,x_CCDF,'-k');
grid on;
title('Tone Reservation');
legend('Orignal','TR');
xlabel('PAPR(dB)');
ylabel('CCDF');
