% 16 QAM 調(diào)制，采用Hadamard變換 N＝128

%OFDM Channel Estimation Based on Comb Pilot
%IFFT_bin_length: IFFT和FFT的點(diǎn)數(shù)
%carrier_count: 子載波個(gè)數(shù)
%bits_per_symbol: 每符號(hào)上的比特?cái)?shù)
%symbols_per_carrier: 每楨的OFDM符號(hào)數(shù)
%Ｘ：欲發(fā)送的二進(jìn)制比特流

clear all;
close all;
N0=128;
IFFT_bin_length0=N0;%ifft長度 prev 128
carrier_count0=N0;%子載波數(shù)   prev 128 
bits_per_symbol=1;%每符號(hào)比特?cái)?shù) prev 2
symbols_per_carrier=100;%一幀符號(hào)數(shù)
N_number0=carrier_count0*symbols_per_carrier*bits_per_symbol;%2400一幀比特?cái)?shù)
M=4;
V=16;
L=4;    %過采樣率
NN=4:.1:12; %CCDF的門限值
GI=8;             % guard interval length
N_snr=100;          % 每比特信噪比
snr=8;             %信噪比間隔

%------------------------------------------------------------
% vector initialization
Xa=zeros(1,N_number0);%2400個(gè)bit

X1=[];
X2=[];
X3=[];
X4=[];
X5=[];
X6=[];
X7=[];
X8=[];
Y1=[];
Y2=[];
Y3=[];
Y4=[];
Y5=[];
Y6=[];
Y7=[];
ccdf0=zeros(1,81);
ccdf1=ccdf0;
ccdf2=ccdf0;
ccdf3=ccdf0;
ccdf4=ccdf0;
ccdf5=ccdf0;
ccdf6=ccdf0;
Xa=randsrc(1,N_number0,[0:3]);%產(chǎn)生二進(jìn)制隨即序列（非0即1）3072

%--------------------------------------------------------
%QAM調(diào)制
X1a=pskmod(Xa,4);

%---------------------------------------------------------
%串并轉(zhuǎn)換
X2a=reshape(X1a,carrier_count0,symbols_per_carrier);%12*128，12個(gè)復(fù)信號(hào)符號(hào)，128列載波
%---------------------------------------------------------
%進(jìn)行空時(shí)分組編碼
N=symbols_per_carrier;
s1=zeros(carrier_count0,symbols_per_carrier);
s2=zeros(carrier_count0,symbols_per_carrier);
s1(:,1:2:symbols_per_carrier)=X2a(:,1:2:symbols_per_carrier);
s1(:,2:2:symbols_per_carrier)=-conj(X2a(:,2:2:symbols_per_carrier));
s2(:,1:2:symbols_per_carrier)=X2a(:,2:2:symbols_per_carrier);
s2(:,2:2:symbols_per_carrier)=conj(X2a(:,1:2:symbols_per_carrier)); 
%---------------------------------------------------------
% oversample L=4
X3a=[s1(1:carrier_count0/2,:);zeros((L-1)*IFFT_bin_length0,symbols_per_carrier);s1(carrier_count0/2+1:carrier_count0,:)];
X3b=[s2(1:carrier_count0/2,:);zeros((L-1)*IFFT_bin_length0,symbols_per_carrier);s2(carrier_count0/2+1:carrier_count0,:)];
%---------------------------------------------------------
% PAPR0 calculation

X4a=ifft(X3a,L*IFFT_bin_length0);  % IFFT Matrix
w0=X4a*carrier_count0 ;
x0=(abs(w0)).^2;
x=x0;
m0=mean(x);
v0=max(x);
papra=10*log10(v0./m0); 

X4b=ifft(X3b,L*IFFT_bin_length0);  % IFFT Transformation
w1=X4b*carrier_count0;
x1=(abs(w1)).^2;
xa=x1;
m1=mean(xa);
v1=max(xa);
paprb=10*log10(v1./m1);
papr0=max(papra,paprb);

%---------------------------------------------------------
% CARI

for k=1:symbols_per_carrier; 
R=randsrc(V,M,[1:4]);   
for v=1:V;
for p=1:M; 
s1a=zeros(N0,1);
s1b=[zeros((p-1)*N0/M,1);-2*s1((p-1)*N0/M+1:p*N0/M,k);zeros((M-p)*N0/M,1)];
s1c=[zeros((p-1)*N0/M,1);-s1((p-1)*N0/M+1:p*N0/M,k);zeros((M-p)*N0/M,1)]+...
  [zeros((p-1)*N0/M,1);s2((p-1)*N0/M+1:p*N0/M,k);zeros((M-p)*N0/M,1)];
s1d=[zeros((p-1)*N0/M,1);-s1((p-1)*N0/M+1:p*N0/M,k);zeros((M-p)*N0/M,1)]-...
  [zeros((p-1)*N0/M,1);s2((p-1)*N0/M+1:p*N0/M,k);zeros((M-p)*N0/M,1)];
s2a=zeros(N0,1);
s2b=[zeros((p-1)*N0/M,1);-2*s2((p-1)*N0/M+1:p*N0/M,k);zeros((M-p)*N0/M,1)];
s2c=[zeros((p-1)*N0/M,1);-s2((p-1)*N0/M+1:p*N0/M,k);zeros((M-p)*N0/M,1)]+...
  [zeros((p-1)*N0/M,1);s1((p-1)*N0/M+1:p*N0/M,k);zeros((M-p)*N0/M,1)];
s2d=[zeros((p-1)*N0/M,1);-s2((p-1)*N0/M+1:p*N0/M,k);zeros((M-p)*N0/M,1)]-...
  [zeros((p-1)*N0/M,1);s1((p-1)*N0/M+1:p*N0/M,k);zeros((M-p)*N0/M,1)];  
if R(v,p)==1
    s1(:,k)=s1(:,k)+s1a;
    s2(:,k)=s2(:,k)+s2a;
elseif R(v,p)==2
    s1(:,k)=s1(:,k)+s1b;
    s2(:,k)=s2(:,k)+s2b;
elseif R(v,p)==3
    s1(:,k)=s1(:,k)+s1c;
    s2(:,k)=s2(:,k)+s2c;
else
    s1(:,k)=s1(:,k)+s1d;
    s2(:,k)=s2(:,k)+s2d;
end
end
X5a(:,k)=ifft(s1(:,k),L*IFFT_bin_length0);
w1a(:,k)=X5a(:,k)*carrier_count0;
x1a(:,k)=(abs(w1a(:,k))).^2;
xaa(:,k)=x1a(:,k);
m1a(:,k)=mean(xaa(:,k));
v1a(:,k)=max(xaa(:,k));
%papr3a(u,k)=10*log10(v1a(:,k)./m1a(:,k));
papr1a(v)=10*log10(v1a(:,k)./m1a(:,k));

X5b(:,k)=ifft(s2(:,k),L*IFFT_bin_length0);
w1b(:,k)=X5b(:,k)*carrier_count0;
x1b(:,k)=(abs(w1b(:,k))).^2;
xab(:,k)=x1b(:,k);
m1b(:,k)=mean(xab(:,k));
v1b(:,k)=max(xab(:,k));
%papr3b(u,k)=10*log10(v1b(:,k)./m1b(:,k));
papr1b(v)=10*log10(v1b(:,k)./m1b(:,k));
papr1(v)=max(papr1a(v),papr1b(v));

end
[papr2(k),Ix(k)]=min(papr1);

%---------------------------------------------------------

%for k=1:symbols_per_carrier;
for l=1:81;
if papr0(k)>NN(l);
    ccdf0(l)=ccdf0(l)+1;
end  

if papr2(k)>NN(l);
    ccdf6(l)=ccdf6(l)+1;
end
end
end

ccdf8=ccdf0./symbols_per_carrier; %N=32
ccdf12=ccdf6./symbols_per_carrier; %N=128,Jacket變換后

NN=4:.1:12;
semilogy(NN,ccdf8,'k',NN,ccdf12,'k--')
title('不同子載波數(shù)的CCDF曲線')
xlabel('papr(dB)'),ylabel('ccdf')
%legend('原來','SLM后',2)
grid on