%-----------------------------------------------------------
% Compute spectrum through nonlinearity
%-----------------------------------------------------------

Nbits   = 10000;
Rs      = 16;     % upsampling rate
rolloff = 0.50;   % rolloff for raised cosine filtering
Ntaps   = 101;    % number of filter taps
Txfilt  = rrcfilt(1/Rs/2,rolloff,201);
BackOff = [1];    % input backoff of nonlinearity

%--- QPSK modulation ----------------
Data = sign(randn(Nbits,1)) + j*sign(randn(Nbits,1));
Data = [1; zeros(Rs-1,1)] * Data';
Data = Data(:);

%--- pi/4-QPSK ----------------------
% Data = sign(randn(1,Nbits)) + j*sign(randn(1,Nbits));
% Rot  = exp(j*pi/4) .^ [1:length(Data)];
% Data = Data .* Rot;
%     
% Data = [1; zeros(Rs-1,1)] * Data;
% Data = Data(:);

%--- OQPSK --------------------------
% Data = [1; zeros(Rs-1,1)] * sign(randn(1,Nbits)) + ... 
%        j * [zeros(Rs/2,1); 1; zeros(Rs/2-1,1)] * sign(randn(1,Nbits));
% Data = Data(:);


TxSig = filter(Txfilt,1,Data); 

Vrms  = std(TxSig);

freq  = Rs*[-512:511]'/1024;  % frequency in symbol rates
temp  = [];

for i=1:length(BackOff)
 	BoSig          = TxSig/Vrms * 10^(-BackOff(i)/20);
	ClipInd        = find(abs(BoSig)>1);
    BoSig(ClipInd) = BoSig(ClipInd) ./ abs(BoSig(ClipInd));
    Nlsig          = BoSig;
    temp           = [temp Nlsig]; 
end

spectra(temp,1024,Rs);