%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% All rights reserved by Krishna Pillai, http://www.dsplog.com 
% The file may not be re-distributed without explicit authorization
% from Krishna Pillai.
% Checked for proper operation with Octave Version 3.0.0
% Author	: Krishna Pillai
% Email		: krishna@dsplog.com
% Version	: 1.0
% Date		: 01 July 2008
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% script for computing the PAPR for a random BPSK/QPSK
% modulation with and without the j rotation

clear all
close all
nFFT = 128; % fft size
nDSC = 114; % number of data subcarriers
nSym = 10^4; % number of symbols

ipBPSK = randsrc(1,nDSC*nSym,[-1 1 ]);  
ipQPSK = 1/sqrt(2)*(randsrc(1,nDSC*nSym,[-1 1 ]) + j*randsrc(1,nDSC*nSym,[-1 1 ]));  
ipModBPSK = reshape(ipBPSK,nDSC,nSym).'; % grouping into multiple symbolsa
ipModQPSK = reshape(ipQPSK,nDSC,nSym).'; % grouping into multiple symbolsa

% Assigning modulated symbols to subcarriers 
xBPSK1F = [zeros(nSym,6) ipModBPSK(:,[1:nDSC/2]) zeros(nSym,3)...
 ipModBPSK(:,[nDSC/2+1:nDSC]) zeros(nSym,5)] ; % without j rotation
xBPSK2F = [zeros(nSym,6) ipModBPSK(:,[1:nDSC/2]) zeros(nSym,3)...
 j*ipModBPSK(:,[nDSC/2+1:nDSC]) zeros(nSym,5)] ; % with j rotation

xQPSK1F = [zeros(nSym,6) ipModQPSK(:,[1:nDSC/2]) zeros(nSym,3)...
 ipModQPSK(:,[nDSC/2+1:nDSC]) zeros(nSym,5)] ;
xQPSK2F = [zeros(nSym,6) ipModQPSK(:,[1:nDSC/2]) zeros(nSym,3)...
 j*ipModQPSK(:,[nDSC/2+1:nDSC]) zeros(nSym,5)] ;

% Taking iFFT, time domain 
xBPSK1t = (nFFT/sqrt(nDSC))*ifft(fftshift(xBPSK1F.')).';
xBPSK2t = (nFFT/sqrt(nDSC))*ifft(fftshift(xBPSK2F.')).';

xQPSK1t = (nFFT/sqrt(nDSC))*ifft(fftshift(xQPSK1F.')).';
xQPSK2t = (nFFT/sqrt(nDSC))*ifft(fftshift(xQPSK2F.')).';

% computing the peak to average power ratio 
meanSquareValueBPSK1 = sum(xBPSK1t.*conj(xBPSK1t),2)/nFFT;
peakValueBPSK1 = max(xBPSK1t.*conj(xBPSK1t),[],2);
paprSymbolBPSK1 = peakValueBPSK1./meanSquareValueBPSK1; 
paprSymbolBPSK1dB = 10*log10(paprSymbolBPSK1);
[nBPSK1 xBPSK1] = hist(paprSymbolBPSK1dB,[0:0.5:15]);

% computing the peak to average power ratio 
meanSquareValueBPSK2 = sum(xBPSK2t.*conj(xBPSK2t),2)/nFFT;
peakValueBPSK2 = max(xBPSK2t.*conj(xBPSK2t),[],2);
paprSymbolBPSK2 = peakValueBPSK2./meanSquareValueBPSK2; 
paprSymbolBPSK2dB = 10*log10(paprSymbolBPSK2);
[nBPSK2 xBPSK2] = hist(paprSymbolBPSK2dB,[0:0.5:15]);

% computing the peak to average power ratio 
meanSquareValueQPSK1 = sum(xQPSK1t.*conj(xQPSK1t),2)/nFFT;
peakValueQPSK1 = max(xQPSK1t.*conj(xQPSK1t),[],2);
paprSymbolQPSK1 = peakValueQPSK1./meanSquareValueQPSK1; 
paprSymbolQPSK1dB = 10*log10(paprSymbolQPSK1);
[nQPSK1 xQPSK1] = hist(paprSymbolQPSK1dB,[0:0.5:15]);

% computing the peak to average power ratio 
meanSquareValueQPSK2 = sum(xQPSK2t.*conj(xQPSK2t),2)/nFFT;
peakValueQPSK2 = max(xQPSK2t.*conj(xQPSK2t),[],2);
paprSymbolQPSK2 = peakValueQPSK2./meanSquareValueQPSK2; 
paprSymbolQPSK2dB = 10*log10(paprSymbolQPSK2);
[nQPSK2 xQPSK2] = hist(paprSymbolQPSK2dB,[0:0.5:15]);

close all
figure
plot(xBPSK1,cumsum(nBPSK1)/nSym,'bs-','LineWidth',2)
hold on
plot(xBPSK2,cumsum(nBPSK2)/nSym,'mp-','LineWidth',2)
plot(xQPSK1,cumsum(nQPSK1)/nSym,'ro-','LineWidth',2)
plot(xQPSK2,cumsum(nQPSK2)/nSym,'g*-','LineWidth',2)
xlabel('papr, x dB')
ylabel('Probability, X <=x')
title('CDF plots of PAPR 128pt FFT (BPSK, QPSK) ')
grid on
legend('BPSK','BPSK-j','QPSK','QPSK-j')
axis([3 12 0 1])