polytime.m

Mtlab toolbox containing useful m files to generate LPI signals

%==========================================================================
% polytime.m
%
%POLYTIME CODE T1-T4
%
clear all;
clc;
disp('*********************************************');
disp('************POLYTIME CODE (T1-T4)************');
disp('*********************************************');

%--------------------------------------------------------------------------
%DEFAULT VARIABLES
%--------------------------------------------------------------------------

A=1;                            % Amplitude of CW
f=1e3;                          % Carrier frequency 
fs=7e3;                         % Sample rate 
SNR_dB = 0;                     % Signal to noise ratio
scale=30;                       % Scaling for plotting time domain graphs
m=2;                            % Number of phase states
SAR=ceil(fs/f);                 % Sampling rate
k=4;                            % Number of stepped frequency segments
T=16e-3;                        % Overall code period
Ts=1/(fs);                      % Sampling period
deltaphi=2*pi/m;                % Fundamental stepsize
deltaf=250;                     % Modulation bandwidth

%--------------------------------------------------------------------------
% NEW INPUT IF NEEDED
%--------------------------------------------------------------------------

newvar = 1;
while newvar == 1;
    disp(' ')
    disp('WHICH PARAMETER DO YOU WANT TO SET ?  ')
    disp(' ')
    fprintf('1. Amplitude of the carrier signal - A= %g.\n', A)
    fprintf('2. Carrier frequency - f (Hz) = %g.\n', f)
    fprintf('3. Sampling frequency - fs (Hz)= %g.\n', fs)
    fprintf('4. Signal to noise ratio - SNR_dB (dB) = %g.\n', SNR_dB)
    fprintf('5. Number of phase states - n = %g.\n', m)
    fprintf('6. Number of segments - k = %g.\n', k)
    fprintf('7. Modulation Bandwidth - df (Hz)= %g.\n', deltaf)
    fprintf('8. Overall Code period - T (s)= %g.\n', T)
    fprintf('9. No changes\n')
    disp(' ')
    option= input('Select a option: ');
    
    switch option
        case 1
            A=input('New amplitude of the carrier signal= '); 
        case 2
            f=input('New carrier frequency (Hz) = ');
        case 3
            fs=input('New sampling frequency (Hz)= ');
        case 4
            SNR_dB=input('New signal to noise ratio (dB)= ');
        case 5
            m=input('New number of phase states =');
        case 6
            k=input('New number of segments =');
        case 7
            deltaf=input('New Modulation Bandwidth =');
        case 8
            T=input('New Overall Code period =');
        case 9
            newvar = 0;
    end
    clc;
end

newvar2 = 1;
while newvar2 == 1;
    disp(' ')
    disp('WHICH PHASECODE DO YOU WANT TO USE ?  ')
    disp(' ')
    disp('1. T1')
    disp('2. T2')
    disp('3. T3')
    disp('4. T4')
    disp(' ')
    option2 = input('Select a phasecode ');
    switch option2
        
        case 1
            
            % Compute the T1 Polytime Phase (formula from paper)
            ttt=1;
            for tt = 0:(Ts):(T-Ts)
                jj = floor(k*tt/T);
                phase(ttt)= mod(((2*pi/m)*floor(((k*tt - jj*T)*(jj*m/T)))), 2*pi);    
                if ttt==1
                    phaseunwrapped(ttt)=phase(ttt);
                else      
                    if phase(ttt)==phase(ttt-1)
                        phaseunwrapped(ttt)=phaseunwrapped(ttt-1);
                    else
                        phaseunwrapped(ttt)=phaseunwrapped(ttt-1)+2*pi/m;
                    end
                end  
                ttt= ttt+1;
            end
            newvar2=0;  
            
        case 2
            
            %Compute the T2 Polytime Phase (formula from paper)
            ttt=1;
            for tt = 0:(Ts):(T-Ts)
                jj = floor(k*tt/T);
                phase(ttt)= mod(((2*pi/m)*floor((((k*tt - jj*T)*((2*jj-k+1)/T)*(m/2))))), 2*pi); 
                if ttt==1
                    phaseunwrapped(ttt)=phase(ttt);
                else
                    if phase(ttt)==phase(ttt-1)
                        phaseunwrapped(ttt)=phaseunwrapped(ttt-1);
                    else
                        phaseunwrapped(ttt)=phaseunwrapped(ttt-1)+2*pi/m;
                    end
                end
                ttt= ttt+1;
            end
            newvar2=0;  
            
        case 3
        
            %Compute the T3 Polytime Phase (formula from paper)
            ttt=1;
            for tt = 0:(Ts):(T-Ts) 
                phase(ttt)=mod(((2*pi/m)*floor((m*deltaf*tt.^2)/(2*T))),2*pi); 
                if ttt==1
                    phaseunwrapped(ttt)=phase(ttt);
                else                      
                    if phase(ttt)==phase(ttt-1)
                        phaseunwrapped(ttt)=phaseunwrapped(ttt-1);
                    else
                        phaseunwrapped(ttt)=phaseunwrapped(ttt-1)+2*pi/m;
                    end
                end  
                ttt = ttt+1;
            end
            newvar2=0;  
            
        case 4
            
            %Compute the T4 Polytime Phase (formula from paper)
            ttt=1;
            for tt = 0:(Ts):(T-Ts)  
                phase(ttt)=mod(((2*pi/m)*floor((m*deltaf*tt.^2)/(2*T)-(m*f*tt)/2)),2*pi);
                if ttt==1
                    phaseunwrapped(ttt)=phase(ttt);
                else 
                    if phase(ttt)==phase(ttt-1)
                        phaseunwrapped(ttt)=phaseunwrapped(ttt-1);
                    else
                        phaseunwrapped(ttt)=phaseunwrapped(ttt-1)+2*pi/m;
                    end
                end 
                ttt = ttt+1;
            end
            newvar2=0;      
    end
    clc;
    
end

Phase=deltaphi*floor(phase/deltaphi);

%--------------------------------------------------------------------------
% This section generates I & Q without T1-T4 phase shift and I & Q with Phase 
% shift.  The signals are generated five times by the outer loop.  The 
% variable 'index' is used to generate a time vector for time domain plots.
%--------------------------------------------------------------------------

index=1;
for p=1:5 %Generate the signal five times and store sequentially in corresponding vectors
    
    for loop=1:(ttt-1) %Loop to shift phase 
        
        I(index)=A*cos(2*pi*f*(loop-1)*Ts+Phase(loop)); %Calculate in phase component of signal with phase shift
        IWO(index)=A*cos(2*pi*f*(loop-1)*Ts); % Calculate in phase component of signal without phase shift
        Q(index)=A*sin(2*pi*f*(loop-1)*Ts+Phase(loop)); % Calculate quadrature component of signal with phase shift
        QWO(index)=A*sin(2*pi*(loop-1)*Ts); %Calculate quadrature component of signal without phase shift
        time(index)=(index-1)*Ts; %time vector cumulation
        index = index+1;
        
    end
    
end

%Power Spectral Density for I with phase shift & with WGN with Signal to noise ratios (SNR) = [0,-5,5,10,-10,-20]
%for loop makes calculations and plots for each value of SNR for WGN
[a,b]=size(I);
SNR=10^(SNR_dB/10);
power=10*log10(A^2/(2*SNR));%calculate SNR in dB for WGN function
noise=wgn(a,b,power);%calculate noise at specified SNR
IN=I+noise;               %add noise to I with P4 phase shift
IPWON=I;                %I with phase shift without noise
QN=Q+noise;            %add noise to Q with P4 phase shift
QPWON=Q;             %Q with phase shift without noise
    

%*******************************************************
%PLOTS
%******************************************************

disp(' ')
plt = input('Do you want to generate plots of the signal (Y/y or N/n) ?','s');
disp(' ')
if (plt == 'Y') | (plt =='y')
    disp(' ')
    
    
    %Plot Power Spectral Density for I without phase shift
    figurecount=1; %figurecount is plot index 
    figure (figurecount); % open new figure for plot
    psd(IWO,[],fs); %Power Spectral Density of I without Phase shift
    title(['Fig # ' num2str(figurecount) ', PSD of I without Phase Shift or Noise']);
    figurecount=figurecount+1; %increment figure count
    
    
    %time domain plot of in phase signal I without phase shift
    figure (figurecount); % open new figure for plot
    % plot small portion of time domain signal I so that data will fit meaningfully in figure.  
    %floor(size(time,2)/scale) selects a small sample of the vectors to plot    
    plot (time,IWO); 
    title(['Fig # ' num2str(figurecount) ', Time Domain of I without Phase Shift or Noise']);
    xlabel('{\itTime - Seconds} '); 
    ylabel('Amplitude');
    grid on;
    figurecount=figurecount+1; %increment figure index
    
    
    %Power Spectral Density for I with phase shift
    figure (figurecount); %open new figure for plot
    psd(I,[],fs); %plot power spectral density of I with phase shift
    title(['Fig # ' num2str(figurecount) ', PSD of I Phase Shift & no Noise']);
    figurecount=figurecount+1; %increment figure index
    %time domain plot of in phase signal I with phase shift
    figure (figurecount); %open new figure for plot
    
     
    % plot small portion of time domain signal I so that data will fit meaningfully in figure.  
    %floor(size(time,2)/scale) selects a small sample of the vectors to plot
    plot (time(1:floor(size(time,2)/scale)),I(1:floor(size(time,2)/scale)));
    title(['Fig # ' num2str(figurecount) ', Time Domain of I with Phase Shift & no Noise']);
    xlabel('{\itTime - Seconds} '); 
    ylabel('Amplitude');
    grid on;
    figurecount=figurecount+1;%increment figure index
    
    
    %Plot PSD and Time Domain of I+ TX Phase + WGN and Time Domain of I + TX Phase    
    figure (figurecount);% open new figure for plot
    psd(IN,[],fs);%plot PSD for specified noise SNR
    title(['Fig # ' num2str(figurecount) ', PSD of I with Phase Shift & Noise SNR=' num2str(10*log10(SNR))]);
    figurecount=figurecount+1;%increment figure index
    
    
    %plot time domain signal I with TX phase shift and WGN at specified SNR
    figure (figurecount);%open new figure for plot
    plot(time(1:floor(size(time,2)/scale)),IN(1:floor(size(time,2)/scale)));
    title(['Fig # ' num2str(figurecount) ', Time Domain of I with Phase Shift & Noise SNR=' num2str(10*log10(SNR))]);
    xlabel('{\itTime - Seconds} '); 
    ylabel('Amplitude');
    grid on;
    figurecount=figurecount+1;%increment figure index
    
    
    figure (figurecount);% open new figure for plot
    plot(time(1:floor(size(time,2)/scale)),IPWON(1:floor(size(time,2)/scale)));
    title(['Fig # ' num2str(figurecount) ', Time Domain of I with Phase Shift witoutout Noise']);
    xlabel('{\itTime - Seconds} ');
    ylabel ('Amplitude');
    grid on;
    figurecount=figurecount+1;%increment figure index
    
    
    % Now check to see if signal is correct by plotting phase shift alone and then determining phase shift from I+jQ.
    % To determine phase shift, look at the phase angle of I+jQ at every 7th time interval.  Expect to see the P4 phase
    % function plot repeated 5 times after unwrapping and detrending the phase angle.
    
    figure(figurecount);%open new figure for plot
    plot(phase);
    title(['Fig # ' num2str(figurecount) ', Phase Shift (radians)']);
    xlabel('i - index for phase change');
    ylabel('Phase Shift - Theta');
    grid on;
    figurecount=figurecount+1;%increment figure index
    %Now strip out points from I and Q to reconstruct phase shift. 
    %I(1:SAR:floor(size(I,2)/5)) selects a data points with the phase values corresponding to the original phase calculation,;
    %by indexing SAR through the first one fifth of the vector computed by floor(size(I,)/5).  The vector was repeated five times.
    signal=I(1:SAR*k:size(I,2))+j*Q(1:SAR*k:size(I,2));
    
    phase_signal=angle(signal);%determine the angle from the complex signal
    
    % unwrap(I) corrects the radian phase angles in array I by adding multiples of