fsk_psk_target.m

Mtlab toolbox containing useful m files to generate LPI signals

% FREQUENCY AND PHASE SHIFT CODE (FSK/PSK)

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

%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
global N;
N=128;                          % Number of fsk/psk sections (frequency hops and phase changes)
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. Sampling frequency - fs (Hz)= %g.\n', fs)
    fprintf('3. Signal to noise ratio - SNR_dB (dB) = %g.\n', SNR_dB)
    fprintf('4. Number of phase/frequency hops - N_C = %g.\n', N)
    fprintf('5. Number of cycles per phase - cpp = %g.\n', cpp)
    fprintf('6. 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
        N=input('New number of phase/frequency hops =');
    case 5
        cpp=input('New number of cycles per phase=');
    case 6
        newvar = 0;           
    end
    clc;
end

airplane; % Compute frequency hops values according to an airplane's frequency response
syn_test; % Generates frequency distribution according to previous probability distribution
    
global target;    % Compute Frequency, first call up airplane variable
global detection;   % 

tb=1/(fs);          % Sampling period
SAR=ceil(fs/max(detection));     % Sampling ratio (for maximum frequency)


% 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; % Time vector for time domain plots.
for ii = (1:N) % Loop for each frequency hop
    f = detection(ii,:); 
        
    for kk = (1:size(target,1))
        if f == target(kk,1)
            init_phase = target(kk,3); %defines initial phase of FSK sequence
        end
    end
    
    phase = rand(1,5*SAR*cpp);    %Compute the phase encoding random sequence for each frequency burst
    greater = find(phase>=0.5);
    phase(greater) = pi;
    lesser = find(phase<0.5);        
    phase(lesser) = 0;                 

    numseq=5;     
    for p=1:numseq %Generate the signal five times and store sequentially in corresponding vectors
        for n=1:SAR*cpp %Loop to increment time for single frequency value.
            I(index+1)=A*cos(2*pi*f*(n-1)*tb+phase(p*n)+init_phase); %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+phase(p*n)+init_phase); % 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


%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(' ')
    global range;
    global airplane_0;
    
    %Plot Target Frequency Distribution and Range/Magnitude plot
    figure;% open new figure for plot
    plot(target(:,1), target(:,2)); title('Ship FFT ABS'); grid on
    title(['Original Target Frequency Probability Distribution']);
    xlabel('Frequency');
    ylabel('Normalized Magnitude = Probability');
    
    figure;% open new figure for plot
    plot(range, airplane_0); title('Ship'); grid on
    title(['Original Target Range / Magnitude Plot']);
    xlabel('Range (ft)');
    ylabel('Magnitude');
    
    % Plot original frequency distribution histogram and frequency random firing distribution
    figure;% open new figure for plot
    orient tall;
    subplot(2,1,1),
    hist(detection(:,1),N);
    %xlabel('Detection Index');
    ylabel('Number of Occurrences');
    fid1=['Target SYNTHETIC'];
    title(fid1);

    subplot(2,1,2),
    bar(target(:,1), target(:,2));
    xlabel('Frequency (Hz)')
    ylabel('Probability')
    fid2=['Target ORIGINAL'];
    title(fid2);
    axis([0 6000 0 0.05]);
    
    %Plot Power Spectral Density for I without phase shift
    figure ; % 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']);
        
    %Power Spectral Density for I with phase shift
    figure ; %open new figure for plot
    psd(I,[],fs); %plot power spectral density of I with phase shift
    title(['PSD of I with Phase Shift & no Noise']);
    
    %time domain plot of in phase signal I with phase shift
    figure ; %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;
        
    %Plot PSD and Time Domain of I+ FSK/PSK Phase + WGN and Time Domain of I + FSK/PSK Phase
    figure ;% 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))]);
    
    %plot time domain signal I with FSK/PSK phase shift and WGN at specified SNR
    figure ;%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;
           
    % 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 last 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_Target_' num2str(ffs) '_' num2str(N) '_' num2str(cpp) '_' num2str(SNR_dB)],'I','Q');
    I=II;
    Q=QQ;
    save(['FSK_PSK_Target_' num2str(ffs) '_' num2str(N) '_' num2str(cpp) '_s'],'I','Q');
    disp(' ');
    disp(['Signal and noise save as :  FSK_PSK_Target_' num2str(ffs) '_' num2str(N) '_' num2str(cpp) '_' num2str(SNR_dB)]);
    disp(['Signal only save as :       FSK_PSK_Target_' num2str(ffs) '_' num2str(N) '_' num2str(cpp) '_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