%  This script shows an example of analyzing a filter/system with the magnitude
%  response, phase response, and group delay.  These are computed directly
%  with basic matlab operations and also using the Matlab's function freqz and
%  grpdelay for comparison.
%
%    Written by Kevin D. Donohue 2/7/2008 ( donohue@engr.uky.edu )
%
fs = 16000;  %  Sampling frequency time domain
df = 1;  % Increment on frequency axis for evaluating spectra
%  Transfer function for analysis 
% 
%              b(1) + b(2)z^-1 + .... + b(m+1)z^-m
%      H(z) = ------------------------------------
%              a(1) + a(2)z^-1 + .... + a(n+1)z^-n
% Coefficient for a near oscillator
b = [0 1];
a = [1.0000   -1.2728    0.8100];
faxis = (0:df:fs/2);  % Frequency axis evaluation points for positive frequencies
z = exp(j*2*pi*faxis/fs);  % Substitute exp(jw) in for z with normalize frequencies 
hn = zeros(size(z));  % Initialze variable to accumulate each term of numerator
%  Accumulate numerator terms
for m=1:length(b)
    hn = hn + b(m)*z.^(-(m-1));
end
hd = zeros(size(z));  % Initialze variable to accumulate each term of denominator
%  Accumulate denominator terms
for n=1:length(a)
    hd = hd + a(n)*z.^(-(n-1));
end
%  Divide numberator and denominator to obtain complex valued transfer
%  function (hd should not be zero for any frequency)
h = hn ./ hd;
%  Plot magnitude response (shows magnitude scaling for any input sinusoid)
figure(1)
plot(faxis,abs(h),'k-','Linewidth', 2)
xlabel('Hz');
ylabel('TF Magnitude')
title('Magnitude Response')
%  Plot phase response (shows phase shifting for any input sinusoid)
figure(2)
plot(faxis,phase(h),'k-','Linewidth', 2)
xlabel('Hz')
ylabel('TF Phase in Radians')
title('Phase Response')
% Compute the group delay
gd = -gradient(phase(h)/(2*pi),df);  % Take gradient of phase values, convert
                              %  phase from radians to Hz (same as frequency increment) 
%  Plot group delay (note one point is lost due to the difference
%  operation)
figure(3)
plot(faxis, gd)
xlabel('Hz')
ylabel('Delay in Seconds')
title('Group Delay')

%   Now apply matlabs functions for the analyses above
[hm, fm] = freqz(b,a,round(fs/(2*df)), fs);  %  Compute Transfer Function
[gdm, fgm] = grpdelay(b,a,round(fs/(2*df)),fs); % Compute Group Delay
figure(4)
plot(fm,abs(hm))
xlabel('Hz');
ylabel('TF Magnitude')
title('Magnitude Response with FREQZ')
figure(5)
plot(fm,phase(hm))
xlabel('Hz')
ylabel('TF Phase in Radians')
title('Phase Response with FREQZ')
figure(6)
%  Note matlab's group delay is in terms of samples not seconds.
plot(fgm,gdm)
xlabel('Hz')
ylabel('Delay in Samples')
title('Group Delay with GRPDELAY')