costas.m

Mtlab toolbox containing useful m files to generate LPI signals

%******************************************************************************
% costas.m
%
% Use:  To generate Costas-coded signals (2 choices)
%
% Inputs:       Amplitude of the signal
%		        Sampling frequency - fs (Hx)
%               Desired SNR (in dB)
%               Costas sequence choice
%
% Output:      In-phase (I) and Quadrature (Q) components of the CW Costas signal
%               (plus plots)
%******************************************************************************

clear all;
clc;
disp('*******************************************');
disp('***************COSTAS CODE ****************');
disp('*******************************************');

%DEFAULT VARIABLES
A=1;                            % Amplitude of signal
fs =15e3;                       % Sampling frequency
SNR_dB = 0;                     % Signal to noise ratio
tp = 1e-3;                      % Sub-period (s)
scale=30;                       % Scaling for plotting time domain graphs
j=sqrt(-1);                     % j

% FREQUENCY CHOICES
newvar = 1;
while newvar == 1;
    disp(' ')
    disp('WHICH FREQUENCY HOP SEQUENCE WOULD YOU LIKE TO USE ?  ')
    disp(' ')
    disp('1.       3, 2, 6, 4, 5, 1 (kHz)');
    disp('2.       5, 4, 6, 2, 3, 1 (kHz)');
    disp('3.       2, 4, 8, 5, 10, 9, 7, 3, 6, 1 (kHz)');
    disp(' ')
    option2= input('Select an option: ');

freq=[3 2 6 4 5 1; 
    5 4 6 2 3 1 ]*1000;

freq_10=[ 2 4 8 5 10 9 7 3 6 1]*1000;

switch option2
    case 1
        seq=freq(1,:);
        [a,length]=size(freq(1,:));
    case 2
        seq=freq(2,:);
        [a,length]=size(freq(2,:));
    case 3
        seq=freq_10(1,:);
        [a,length]=size(freq_10(1,:));
end
    newvar=0;
    clc;
end

% NEW INPUT 
newvar = 1;
while newvar == 1;
    disp(' ')
    disp('WHICH COSTAS CW FREQUENCY HOPPING SIGNAL PARAMETER DO YOU WANT TO SET ?  ')
    disp(' ')
    fprintf('1. Amplitude - A= %g.\n', A)
    fprintf('2. Sampling frequency - fs (Hz)= %g.\n', fs)
    fprintf('3. Signal to noise ratio - SNR_dB (dB) = %g.\n', SNR_dB)
    fprintf('4. Frequency duration (s)  = %g.\n', tp)
    fprintf('5. No changes\n')
    disp(' ')
    option= input('Select a option: ');
    
    switch option
    case 1
        A=input('New amplitude of the signal= '); 
    case 2
        fs=input('New sampling frequency (Hz) ');
    case 3
       SNR_dB=input('New signal to noise ratio (dB)= ');
    case 4
       tp=input('New duration of each frequency = ');
    case 5
        newvar = 0;
    end
    clc;
end
    
tb=1/(fs);                   % Sampling period

% This section generates I & Q without COSTAS 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=0;
numseq=5;
for p=1:numseq %Generate the signal five times and store sequentially in corresponding vectors
    for xx=1:length
        cpf=tp*seq(xx);         %number of cycles of frequency within the code period tp
        SAR=ceil(fs/seq(xx));   %number of samples needed per cycle
            for n=1:SAR*cpf
            I(index+1)=A*cos(2*pi*seq(xx)*(n-1)*tb); 
            Q(index+1)=A*sin(2*pi*seq(xx)*(n-1)*tb); 
            time(index+1)=index*tb; %time vector cumulation
            index = index +1;
            end
    end
end        


%Power Spectral Density for I with phase shift & with WGN with Signal to noise ratios 
%for loop makes calculations and plots for each value of SNR for WGN
[a,b]=size(I);
samps_seq=b/numseq; %Samples in a sequence

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 COSTAS phase shift
IPWON=I;                %I with phase shift without noise
QN=Q+noise;            %add noise to Q with COSTAS 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(' ')
    
    figurecount=1;
    %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

    
    
    %Plot PSD and Time Domain of I+ COSTAS Phase + WGN and Time Domain of I 
    
    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
    
   
   
else
    disp('Signal not plotted')
    fprintf('\n\n')
end
    
% This section generates the files for analysis

INP=IN';%transpose I with noise and COSTAS phase shift for text file
QNP=QN';%transpose Q with noise and COSTAS phase shift for text file
IPWONT=IPWON';%transpose I with phase without noise for text file
QPWONT=QPWON';%transpose Q with phase without noise for text file

% % save results in data files

I= INP(:,1);
Q=QNP(:,1);
 
II= IPWONT(:,1);
QQ=QPWONT(:,1);
 
disp(' ')
saveresult = input('Do you want to save the new signal (Y/y or N/n) ?','s');

if (saveresult == 'Y') | (saveresult =='y') 
    ff=7;
    ffs=floor(fs/1e3);
    tpp=tp/1e-6;
    save(['C_' num2str(option2) '_' num2str(ffs) '_' num2str(tpp) '_' num2str(SNR_dB)],'I','Q');
    I=II;
    Q=QQ;
    save(['C_' num2str(option2) '_' num2str(ffs) '_' num2str(tpp) '_s'],'I','Q');
    disp(' ');
    disp(['Signal and noise save as :  C_' num2str(option2) '_' num2str(ffs) '_' num2str(tpp) '_' num2str(SNR_dB)]);
    disp(['Signal only save as :         C_' num2str(option2) '_' num2str(ffs) '_' num2str(tpp) '_s']);
    disp(['Directory:                       ' num2str(cd)]); 
    disp(['NOTE: Number of sequences = ' num2str(numseq) ' Samples in single FH sequence = ' num2str(samps_seq)]);
else
    disp(' ')
    disp('Signal not saved')
    fprintf('\n\n')
end