%==========================================================================
% The mfile investigates the effects of high power amplifier on the ofdam
% signals. The effects on spectrum ang Modulation Error Rate (MER) is of
% more concern.
% 
% Written By     : Hamid Ramezani
% Date           : 17-Jun-2007
% Code version   : 1
% Matlab Version : 7.4.0.287 (R2007a)
%==========================================================================

% Initialization
    clear all;
    close all;
    clc;
    
%==========================================================================
% Setting Parameters
%==========================================================================
    % OFDM System Parameters
    N           = 256;      % length of OFDM IFFT (16,32,64,...,2^n)
    M           = 64;       % number of QAM constellation points (4,16,64,256)
    numOfZeros  = N/4+1;    % numOfZeros must be an odd number and lower 
                            % than N. The zero padding operation is
                            % necessarry in practical implementations.
    GI          = 1/4;      % Guard Interval (1/4,1/8,1/16,...,4/N)
    BW          = 20;       % OFDM signal Band width in MHz
    numOfSym    = 100;      % number of OFDM Symbols
    
    % Amplifire Parameters
    satLevel = 5;   % in dB , higher than the tx out mean of voltage
    
%==========================================================================
%   Main Program
%==========================================================================
    txData     = randint(N-numOfZeros,numOfSym,M);  % data generation
    
    % QAM modulation
    txDataMod  = qammod(txData,M);
    
    % zeros padding
    txDataZpad = [txDataMod((N-numOfZeros+1)/2:end,:);...
                  zeros(numOfZeros,numOfSym);...
                  txDataMod(1:(N-numOfZeros+1)/2+1,:)];
        % Note : in practice zero padding operation must be followed by
        % a standard. Usually the last part of data frame shifts to the first
        % part of zero padded frame.
        
   % IFFT
    txDataZpadIfft = sqrt(N)*ifft(txDataZpad,N);
   
    % Guard Interval Insertion
    txDataZpadIfftGI    = [txDataZpadIfft((1-GI)*N+1:end,:);txDataZpadIfft];
    
    % Amplifier Model
    txDataZpadIfftGIAbs     = abs(txDataZpadIfftGI);            % tx data amplitude

    % tx data amplitude standard deviation and mean
    txDataZpadIfftGIAbsStd  = mean(std(txDataZpadIfftGIAbs));   
    txDataZpadIfftGIAbsMean = mean(mean(txDataZpadIfftGIAbs));
    
    % tx data phase in radian 
    txDataZpadIfftGIAng     = angle(txDataZpadIfftGI);
    
    % It is imagined that the amplifier has no effect on the phase of the
    % signal. The solid state amplifier effects on signal phase is about 5
    % degrees. The amplifier AM/AM response is followed by x/sqrt(1+(x/k)^2)
    txDataZpadIfftGIAbsHPA = txDataZpadIfftGIAbs ./...
           sqrt(1+(txDataZpadIfftGIAbs/(txDataZpadIfftGIAbsMean*10^(satLevel/10))).^2);
    % no change in the phase
    txDataZpadIfftGIAngHPA = txDataZpadIfftGIAng;
    
    % mean of amplitude after amplification
    txDataZpadIfftGIAbsHPAmean = mean(mean(txDataZpadIfftGIAbsHPA));
    % standard deviation after amplification
    txDataZpadIfftGIAbsHPAStd  = mean(std(txDataZpadIfftGIAbsHPA));
    

    % polar to cartesian conversion
    txDataZpadIfftGIHPA    = txDataZpadIfftGIAbsHPA.* ...
                            exp(sqrt(-1) * txDataZpadIfftGIAngHPA);                            

    % receiver part
    % Guard Interval removal
    rxDataZpadIfftHPA  = txDataZpadIfftGIHPA(GI*N+1 : N+GI*N,:);
    % FFT operation
    rxDataZpadHPA      = 1/sqrt(N)*fft(rxDataZpadIfftHPA,N);
    % zero removal and rearrangement
    rxDataModHPA       = [rxDataZpadHPA((N-(N-numOfZeros-1)/2+1):N,:);...
                          rxDataZpadHPA(1:(N-numOfZeros+1)/2,:)];
    % demodulation
    rxDataHPA          = qamdemod(rxDataModHPA/mean(std(rxDataModHPA))*mean(std(txDataMod)),M);
    
%==========================================================================
%       statistical computation
%==========================================================================
    % Mean Error Rate computation
    MER       = 10*log10(mean(var(rxDataModHPA./mean(std(rxDataModHPA))...
                 - txDataMod./mean(std(txDataMod)))));
    % Bit Error Rate computation
    [num BER] = symerr(rxDataHPA,txData);
%==========================================================================
% graphical observation
%==========================================================================
   f1 = figure(1);    
   set(f1,'color',[1 1 1]);
        subplot(2,2,1);
            
            spectrumFftSize = 2*N;
            % spectrum of signal befor High Power Amplifier
            txSpec  = 20*log10(mean(abs(fft(txDataZpadIfftGI(:,:)./ ...
                      mean(std(txDataZpadIfftGI)),spectrumFftSize)),2));
            % spectrum of signal after High Power Amplifier
            HpaSpec = 20*log10(mean(abs(fft(txDataZpadIfftGIHPA(:,:)./ ...
                      mean(std(txDataZpadIfftGIHPA)),spectrumFftSize)),2));
            % corresponding frequency 
            Freq    = linspace(-BW/2,BW/2,length(txSpec));
            
            plot(Freq,[txSpec(length(txSpec)/2:length(txSpec));...
                txSpec(1:(length(txSpec)/2-1))]);
            hold on
            grid on
            plot(Freq,[HpaSpec(length(txSpec)/2:length(txSpec));...
                HpaSpec(1:(length(txSpec)/2-1))],'r');
            grid on;
            xlabel('representing frequency');
            ylabel('spectrum signals (first symbol)');
            title('Spectrum Effects')
            legend('Befor Amplifier','After Amplifier')

        subplot(2,2,2);
            plot(real(reshape(rxDataModHPA,1,numOfSym*(N-numOfZeros)))/mean(std(rxDataModHPA)),...
                 imag(reshape(rxDataModHPA,1,numOfSym*(N-numOfZeros)))/mean(std(rxDataModHPA)),'.r')
            hold on
            plot(real(reshape(txDataMod,1,numOfSym*(N-numOfZeros)))/mean(std(txDataMod)),...
                 imag(reshape(txDataMod,1,numOfSym*(N-numOfZeros)))/mean(std(txDataMod)),'.b')
            xlabel('I channel');
            ylabel('Q channel');
            title('signal Constelleations');
            legend('After Amplifier','Before Amplifier');

    	subplot(2,2,3)
            % normalize amplitude 
            txAmp = linspace(0,5,100);
            amAmp = txAmp./(1+(txAmp/satLevel).^2);
            plot(txAmp,txAmp,'b');
            hold on
            plot(txAmp,amAmp,'r');
            plot(1,1,'om');
            xlabel('Input Amplitude');
            ylabel('Output Amplitude');
            title('AM/AM response of power amplifier');
            legend('Linear Response','Amplifier Response','Mean of OFDM Amplitude');
            
        subplotHandel = subplot(2,2,4);
            text(0,1 ,['Mean Error Rate     : ',num2str(MER),' dB']);
            text(0,.8,['Bit  Error Rate     : ',num2str(BER)]);
            
            text(0,.6,['Modulation          : ',num2str(M),' QAM']);   
            text(0,.4,['IFFT Size           : ',num2str(N),' points']);  
            text(0,.2,['Guard Interval Size : ',num2str(N),' points']);              
            text(0,.0,['Saturation Level    : ',num2str(satLevel),' dB relative to AM Avg']);
            % setting the axes invisibale
            set(subplotHandel,'Xcolor',[1 1 1]);
            set(subplotHandel,'Ycolor',[1 1 1]);