% Compute energy spectra for Impulse, Doublet, 3-2-1-1 and Modified 3-2-1-1
%
% This program computes power spectral densities of input signals
% using Eq. (2.8) of Chapter 2:
%
% Using the function "EnergySpectrum" a spectrum of any general multistep
% input can be calculated and plotted.
% The following input parameters are necessary for EnergySpectrum:
% Inputs: 
%     V        Vector with amplitudes of step signals (both, amplitudes and number
%              of steps are defined by V)
%     dtSig    Length of one time step in seconds
%     wa       Starting frequency
%     we       Final frequency
%     nw       Number of frequency points (default=100)
%
% Outputs:
%     w        Frequencies [1/s]
%     E        Energy spectrum
% 
% Chapter 2: Data Gathering
% "Flight Vehicle System Identification - A Time Domain Methodology"
% Author: Ravindra V. Jategaonkar
% published by AIAA, Reston, VA 20191, USA

clear all;
close all;

%--------------------------------------------------------------------------
% 3211-Signal
dt3211 = 0.3;                                      % time step (1 sec dt)
V      = [0.3 0.3 0.3 -0.3 -0.3 0.3 -0.3];         % amplitudes
[w3,E3] = EnergySpectrum(V,dt3211,0.01,20,100);

%--------------------------------------------------------------------------
% Modified 3211-signal
dtM3211 = 0.3;                                     % time step
V       = [0.24 0.24 0.24 -0.36 -0.36 0.33 -0.33]; % amplitudes (0.8, 1.2, 1.1, 1.1)
[w4,E4] = EnergySpectrum(V,dtM3211,0.01,20,100);

%--------------------------------------------------------------------------
% Doublet
% The time step of doublet is obtained as follows:
% Assuming a dt3211 = 0.3, leads to omega*dt3211 = 1.6 (Eq. 2.11), which means
% omega = 1.6/0.3. Now from Eq. (2.9), Omega*dt11 = 2.3, the time step for the
% doublet input is given by: 2.3*0.3/1.6 = 0.4312.
dt11 = 0.43;                                       % time step
V    = [0.6 -0.6];                                 % amplitudes
[w2,E2] = EnergySpectrum(V,dt11,0.01,20,100);

%--------------------------------------------------------------------------
% Pulse
dt1 = 2.1;                                         % time step (same total duration)
V   = [1.0];                                       % amplitude
[w1,E1] = EnergySpectrum(V,dt1,0.0001,20,100);

%--------------------------------------------------------------------------
% plot
wn = w3*dt3211;
plot(wn,E1/dt1^2,'b', wn,E2/dt11^2,'r', wn,E3/dt3211^2,'k', wn,E4/dtM3211^2,'m'); grid;
axis([0 5, 0 1.2]);
xlabel('normalized frequency dtsmall'); ylabel('E/\Deltat2');
legend('Impulse', 'Doublet', '3211', 'Modified 3211', 1);