fsk_psk_costas.m

Mtlab toolbox containing useful m files to generate LPI signals

% FREQUENCY AND PHASE SHIFT CODE using Costas Frequency modulation (FSK/PSK)

clear all;
clc;
disp('*******************************************************************************************');
disp('*************FREQUENCY AND PHASE SHIFT CODE (FSK/PSK USING COSTAS/BARKER CODING *************');
disp('*******************************************************************************************');

tp=1e-3; % Frequency duration

%DEFAULT VARIABLES
A=1;                          % Amplitude of CW
fs =15000;                    % Sampling Frequency 
SNR_dB = 0;                   % Signal to Noise Ratio
scale=100;                    % Scaling for plotting time domain graphs
j=sqrt(-1);                    % j
barker = 5;                   % Number of bits for Barker Code for phase modulation  

% 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. Sampling frequency - fs (Hz)= %g.\n', fs)
    fprintf('3. Signal to noise ratio - SNR_dB (dB) = %g.\n', SNR_dB)
    fprintf('4. Number of bits per Barker code for phase modulation - barker (13/11/7/5)= %g.\n', barker)
    fprintf('5. No changes\n')
    disp(' ')
    option= input('Select a option: ');
    
    switch option
    case 1
        A=input('New amplitude of the carrier signal= '); 
    case 2
        fs=input('New sampling frequency (Hz)= ');
    case 3
        SNR_dB=input('New signal to noise ratio (dB)= ');
    case 4
        barker=input('New number of bits for Barker Code = ');  
    case 5
        newvar = 0;           
    end
    clc;
end

% FREQUENCY CHOICES
newvar = 1;
while newvar == 1;
    disp(' ')
    disp('WHICH FREQUENCY 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(' ')
    option2= input('Select an option: ');

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

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

%maximum=max(seq);
tb=1/fs;

%cpp=11;
% This section generates I & Q without phase shift and I & Q with Phase shift.  The signals are generated
% five times the number of frequency hops 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;
numseq=2;
for p=1:numseq %Generate the signal five times and store sequentially in corresponding vectors
    
        for xx=1:length
        
        cpf=tp*seq(xx);         %cycles per frequency within the period tp
        SAR=ceil(fs/seq(xx));   %
        
 % Create spanning Barker frequency
clear phase;
if barker==13
    phase = [zeros(1,SAR*cpf*5),(ones(1,SAR*cpf*2)*pi),zeros(1,SAR*cpf*2),(ones(1,SAR*cpf)*pi),zeros(1,SAR*cpf),(ones(1,SAR*cpf)*pi),zeros(1,SAR*cpf)];% 13 bits
elseif barker==11
    phase = [zeros(1,SAR*cpf*3),(ones(1,SAR*cpf*3)*pi),zeros(1,SAR*cpf),(ones(1,SAR*cpf)*pi),(ones(1,SAR*cpf)*pi),zeros(1,SAR*cpf),(ones(1,SAR*cpf)*pi)];% 11 bits
elseif barker==7
    phase = [zeros(1,SAR*cpf*3),(ones(1,SAR*cpf*2)*pi),zeros(1,SAR*cpf),ones(1,SAR*cpf)*pi];% 7 bits
elseif barker==5
    phase = [zeros(1,SAR*cpf*3),(ones(1,SAR*cpf)*pi),zeros(1,SAR*cpf)];% 5 bits
end
%phase=phase.*pi;
             
            for n=1:SAR*cpf*barker
        
                IWO(index+1)=A*cos(2*pi*seq(xx)*(n-1)*tb); % Calculate in phase component of signal without phase shift
    
                QWO(index+1)=A*sin(2*pi*seq(xx)*(n-1)*tb); %Calculate quadrature component of signal without phase shift
                
                I(index+1)=A*cos(2*pi*seq(xx)*(n-1)*tb + phase(n) - pi/2.2); %Calculate in phase component of signal with phase shift
                
                Q(index+1)=A*sin(2*pi*seq(xx)*(n-1)*tb + phase(n) - pi/2.2); % Calculate quadrature component of signal with phase shift
    
                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 (SNR) = [0,-5,5,10,-10,-20]
%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 FSK/PSK phase shift
IPWON=I;                %I with phase shift without noise
QN=Q+noise;            %add noise to Q with FSK/PSK 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
     figure ; % open new figure for plot
     psd(IWO,2048,fs); %Power Spectral Density of I without Phase shift
     title(['PSD of I without Phase Shift, without Noise']);
     axis([0 fs/2 -30 30]);

    %Power Spectral Density for I with phase shift
    figure ; %open new figure for plot
    psd(IPWON,2048,fs); %plot power spectral density of I with phase shift
    title(['PSD of I with Phase Shift & No Noise']);
        axis([0 fs/2 -30 30]);
    %Power Spectral Density for I with phase shift
    figure ; %open new figure for plot
    psd(IN,2048,fs); %plot power spectral density of I with phase shift
    title(['PSD of I with Phase Shift & Noise SNR=' num2str(10*log10(SNR))]);
         axis([0 fs/2 -30 30]);
           
    % 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 FSK/PSK phase
    % function plot repeated 5 times after unwrapping and detrending the phase angle.
    
    figure;%open new figure for plot
    plot(phase);
    title(['FSK/PSK Phase Shift for each frequency hop']);
    xlabel('i - index for phase change');
    ylabel('FSK/PSK Phase Shift - Theta');
    grid on;
else
    disp('Signal not plotted')
    fprintf('\n\n')
end
    
% This section generates the files for analysis

INP=IN';%transpose I with noise and FSK/PSK phase shift for text file
QNP=QN';%transpose Q with noise and FSK/PSK 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') 
    ffs=floor(fs/1e3);
    save(['FSK_PSK_Costas_' num2str(ffs) '_' num2str(barker) '_' num2str(SNR_dB)],'I','Q');
    I=II;
    Q=QQ;
    save(['FSK_PSK_Costas_' num2str(ffs) '_' num2str(barker) '_s'],'I','Q');
    disp(' ');
    disp(['Signal and noise save as :  FSK_PSK_Costas_' num2str(ffs) '_' num2str(barker) '_' num2str(SNR_dB)]);
    disp(['Signal only save as :       FSK_PSK_Costas_' num2str(ffs) '_' num2str(barker) '_s']);
    disp(['Directory:                       ' num2str(cd)]); 
    disp(['NOTE: Number of sequences = ' num2str(numseq) ' Samples in single sequence = ' num2str(samps_seq)]);
else
    disp(' ')
    disp('Signal not saved')
    fprintf('\n\n')
end