% Quick way to make bode plots and a bunch of other useful plots to fully
% characterize filter performance.
function bodeplot(f, G, color)
% f = an array of frequencies.
% G = the corresponding complex response.
% color = optional specification for the appearance of the line, e.g. '-r'.

gaindB = 20*log10(abs(G));
phase = angle(G)*180/pi; %in degrees

%%%%%%%%%%%%%%%%%%%%%%%%%
%  B O D E   P L O T S  %
%%%%%%%%%%%%%%%%%%%%%%%%%
% Gain plot
figure(1)
semilogx(f, gaindB, color, 'linewidth', 6)
hold on

% Phase plot
figure(2)
semilogx(f, phase, color, 'linewidth', 6)
hold on

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  T I M E   D E L A Y   P L O T S  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
timedelay = angle(G)./(2*pi*f);

% Time delay plot
figure(3)
semilogx(f, timedelay, color, 'linewidth', 6)
hold on

% Time delay vs. attenuation
figure(4)
plot(-gaindB, timedelay, color, 'linewidth', 6)
hold on

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  S T E P   &   I M P U L S E   R E S P O N S E  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Here we assume that there is enough information in the tabulated response
% function, G(f), to work out a step response. But this must involve the FFT,
% with fixed stepsize in frequency (rather than logarithmic spacing, as
% would be more typical for a Bode plot). So we interpolate onto a frequency
% axis with fixed step size.

% Time axis
N = 4096; % Best to choose a power of two, to speed up the FFT. Also need
   % a number large enough to resolve all the detail. This can be particularly
   % challenging with the impulse response.
time = linspace(-5, 15, N); % Assume that 1s is the relevant timescale;
   % This allows 10 time units in each half of the square wave defined below.
   % The actual step will run from 0-10 seconds.
dt = time(2)-time(1); % Interval on time axis

% Step function
step(1 : N/4)  = 0;
step(N/4+1 : 3*N/4) = 1;
step(3*N/4+1 : N) = 0;
fstep = fft(step); 

% Impulse function
impulse(1:N) = 0;
impulse(N/4) = 1/dt; % A properly normalized delta function
fimpulse = fft(impulse);

% Frequency axis ff (arranged as it should be for an FFT)
ff1 = (0:(N/2))/(N*dt); % Positive frequencies, of which there are N/2 + 1.
ff2 = -( (N/2 - 1) : -1 : 1)/(N*dt); % Negative frequencies, "  "  N/1 - 1.
ff = [ff1, ff2]; 
   % The full frequency array never actually gets used, but here it is.

% Response function GG for all values of the frequency ff
GG1 = interp1(f, G, ff1, 'spline', 'extrap'); % for the positive frequencies.
GG2 = interp1(f, conj(G), -ff2, 'spline', 'extrap'); % for the negative frequencies.
   % Since usually only positive frequencies have been tabulated in argument f,
   % I've taken advantage of the fact that G is known to be the Fourier
   % transform of a real function.
GG = [GG1, GG2]; % Response arranged like a proper FFT.

% Calculate step response
stepresponse = real( ifft(GG.*fstep) );
impulseresponse = real( ifft(GG.*fimpulse) );

% Plotting step and impulse responses
figure(5)
plot(time, stepresponse, color, 'linewidth', 6)
hold on

figure(6)
plot(time, impulseresponse, color, 'linewidth', 6)
hold on
