%  This script illustrates plotting the frequency response given a H(z)
%  where  H(z) = 1/(1-.5*z^-1) directly and using Matlab's FREZQ command.
%  The impuls resoponse is also computed and plotted using Matlab's 
%  Filter Command.
%
%   Updated by Kevin D. Donohue (donohue@engr.uky.edu) January 2010

clear
npts = 256;  % number of points in plot
% f = (-fix(npts/2):(fix(npts)/2)-1)/npts;  %  Create normalized frequency axis (double sided)
f = .5*(0:npts-1)/npts;  %  Create normalized frequency axis (positive only)
w = 2*pi*f;  %  Convert Hz to radians per second

 z = exp(1j*w); %  Create z vector of samples along the unit circle
f0 = 1/8;   z0 = exp(-1j*2*pi*f0); %  Excitation frequency in Hz 

%  Evaluate points along the unit circle of this H(z)
hz = 1 ./ (1 - .5*z.^-1);  
figure(1); plot(f,abs(hz))  %  Plot magnitude
hz0 = 1 ./ (1 - .5*z0^-1); % Evaluate TF magnitude point at f0
ahz = abs(hz0);
title(['Magnitude of H(z), at f = ' num2str(f0) ' Hz the magnitude is ' num2str(ahz)])
xlabel('Normalize Hertz')
ylabel('Amplitude')
figure(2); plot(f,180*phase(hz)/pi)  % Plot Phase
phz = angle(hz0);
title(['Phase of H(z), at f = ' num2str(f0) ' Hz, the phase is ' num2str(phz) ' radians'])
xlabel('Normalize Hertz')
ylabel('Degrees')

%   Plot impulse response
%
ximp = [1, zeros(1,16)];  %  Create impulse input
%  Get filter coefficient from Hz corresponding to the difference equaiton
a = [1 -.5];   %  Filter coefficients (denominator of hz)
b = 1;  %  Filter coefficients (numerator of hz)
yimp = filter(b,a,ximp);
figure(3); plot((0:length(yimp)-1),yimp)
xlabel('Samples')
ylabel('Amplitudes')
title('Impulse Response')

% Plotting Frequency response using matlabs FREQZ function
fs= 1;  % Normalize Sampling Frequency
[hzz, fax] = freqz(b, a, npts, fs);
figure(4); plot(fax,abs(hzz))  %  Plot magnitude 
title('Magnitude of H(z) with freqz function')
xlabel('Normalize Hertz')
ylabel('Amplitude')
figure(5); plot(fax,180*angle(hzz)/pi)  % Plot Phase
title('Phase of H(z) with freqz function')
xlabel('Normalize Hertz')
ylabel('Degrees')