frank.m

Mtlab toolbox containing useful m files to generate LPI signals

%******************************************************************************
% frank.m
%
% Use:  Code to generated FRANK-coded signals
%
% Inputs:       Amplitud of the carrier signal
%		        Carrier signal frequency - f (Hz)
%               Sampling frequency - fs (HZ)
%               Desired SNR in dB
%               Number of phases
%               Number of cycles per phase
%               
%
%Output:        In-phase (I) and Quadarture (Q) components of the signal
%               Plots
% To change number of code periods generated (currently 5) 
% go to line 89 and change p=1:5
%******************************************************************************

clear all;
clc;
disp('*******************************************');
disp('**************FRANK  CODE ******************');
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
j=sqrt(-1);                   % j
m=8;                          % Number of code phases
cpp = 1;                      %Number of cycles per phase

% NEW INPUT 
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 codes - m = %g.\n', m)
    fprintf('6. Number of cycles per phase - cpp = %g.\n', cpp)
    fprintf('7. 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 codes =');
    case 6
        cpp=input('New number of cycles per phase=');
    case 7
        newvar = 0;
    end
    clc;
end

SAR=ceil(fs/f);          % Sampling ratio
tb=1/(fs);                   % Sampling period


N=m;
%Creating the phase matrix

    for i=1:m
        for j=1:m
            phi(i,j)=2*pi/N*(i-1)*(j-1);
        end
    end
    
% This section generates I & Q without phase shift and I & Q with Frank 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. 
% The signal is generated at seven samples per phase change.   

index=0;
for p=1:5 %Generate the signal five times and store sequentially in corresponding vectors
    for i=1:m %Loop to shift phase 
        for j=1:m
            for n=1:SAR*cpp                                           %Loop to increment time for single phase value.
                I(index+1)=A*cos(2*pi*f*(n-1)*tb+phi(i,j)); %Calculate in phase component of signal with phase shift
                IWO(index+1)=A*cos(2*pi*f*(n-1)*tb); % Calculate in phase component of signal without phase shift
                Q(index+1)=A*sin(2*pi*f*(n-1)*tb+phi(i,j)); % Calculate quadrature component of signal with phase shift
                QWO(index+1)=A*sin(2*pi*f*(n-1)*tb); %Calculate quadrature component of signal without phase shift
                time(index+1)=index*tb; %time vector cumulation
                index = index +1;
            end
        end
    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 Frank phase shift
IPWON=I;                %I with phase shift without noise
QN=Q+noise;            %add noise to Q with Frank phase shift
QPWON=Q;             %Q with phase shift without noise
    
ff=floor(f/1e3);
ffs=floor(fs/1e3);
%*******************************************************
%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([' 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(1:floor(size(time,2)/scale)),IWO(1:floor(size(time,2)/scale))); 
    title([ ' 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([' 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([' 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+ Frank Phase + WGN and Time Domain of I + Frank Phase
    
    figure (figurecount);% open new figure for plot
    psd(IN,[],fs);%plot PSD for specified noise SNR
    title([' PSD of I with Phase Shift & Noise SNR=' num2str(10*log10(SNR))]);
    figurecount=figurecount+1;%increment figure index
    
    %plot time domain signal I with Frank 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([' 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([' 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 Frank phase
    % function plot repeated 5 times after unwrapping and detrending the phase angle.
    
    figure(figurecount);%open new figure for plot
    plot(phi);
    title(['Frank Code Phase Shift']);
    xlabel('i - index for phase change');
    ylabel('Frank Phase Shift - Theta');
    grid on;
    figurecount=figurecount+1;%increment figure index
    
    
    figure(figurecount);%open new figure for plot
    
   %to concatenate the phase matrix in only 1-row matrix
    nn=0;
    for ii=1:N
        for jj=1:N
            nn=nn+1;
            phi2(nn)=phi(ii,jj);
        end
    end
    
    xx=0:length(phi2)-1;
    stairs(xx,phi2);grid
    title(['']);
    xlabel('i - index for phase change');
    ylabel('Frank phase shift (rad)');
    figurecount=figurecount+1;%increment figure index
    
    figure(figurecount);%open new figure for plot
    undoo=rem(phi2, 2*pi); 
    stairs(undoo); grid;
    ylabel('Signal phase (rad)');
    xlabel('i - index for phase change');
    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*cpp:size(I,2))+j*Q(1:SAR*cpp: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