%FILTDEMO Demonstration of filter design techniques.
echo on
clc
%	This example demonstrates the use of the YULEWALK, BUTTER and
%	CHEBY1 functions to design bandpass filters.  A sampling rate
%	of 1000 Hz is assumed.

%	The YULEWALK function allows you to specify a piecewise shape
%	for the desired frequency response magnitude.  YULEWALK then
%	finds an infinite-impulse response filter of the desired order
%	that fits the frequency response in a least-squares sense.
%
%	We start by specifying the desired frequency response, point-wise,
%	with 1.0 corresponding to half the sample rate:

f = [0 .4 .4 .6 .6 1];
H = [0  0   1   1   0  0];

pause % Strike any key to continue.

clc
%	Next we plot the desired frequency response to make sure it is what
%	we want.  To plot the frequency response, we can scale the
%	normalized frequency axis:

fs = 1000;		% sampling rate
fhz = f*fs/2;
  
%	Now we are ready to plot.

pause % Strike any key for plot.

plot(fhz,H), title('Desired Frequency Response'), ...
xlabel('Frequency (Hz)'), ylabel('Magnitude'),  pause
clc
% 	If the desired frequency response did not look right, 
% 	we would go back and change H. However, in this case
% 	everything looks OK so we perform the design:

N = 8;		% Order of the filter (number of poles and zeros).

[Bh,Ah] = yulewalk(N,f,H);		% Working, please wait.....

%	The design is complete and we can look at the filter coefficients.

pause % Strike any key to continue.
format compact
clc
Ah,Bh	% Here are the filter coefficients:
pause % Strike any key to continue.
format
clc
%	Of course it is impossible to know by looking at the coefficients
%	whether the design is good or bad. Let's compute the frequency
%	response of the filter:

n = 256;
hh = freqz(Bh,Ah,n);	% compute complex frequency response
hy  = abs(hh);    	% compute magnitude

%	Now we can plot the frequency response magnitude
%	and compare it to the desired response:

pause % Strike any key for plot.
ff  = fs/(2*n) * (0:n-1);
plot(fhz,H,ff,hy), title('Actual vs. Desired Frequency Response'), ...
xlabel('Frequency (Hz)'), ylabel('Magnitude'), pause
clc
%	Now let's design Butterworth and Chebyshev bandpass filters
%	with the same passband (defined between 0.0 and 1.0).

N = 4;				% Filter order
passband = [.4 .6];		% Passband specification
[Bb,Ab] = butter(N, passband);

ripple = .1;			% Allowable ripple, in decibels
[Bc,Ac] = cheby1(N, ripple, passband);

% 	Now find the frequency responses:

hb = freqz(Bb,Ab,n);
hc = freqz(Bc,Ac,n);
h = [abs(hh) abs(hb) abs(hc)];

pause % Strike any key for plot.

plot(ff,h), title('YuleWalk, Butterworth and Chebyshev filters'), ...
xlabel('Frequency (Hz)'), ylabel('Magnitude'),pause
clc
%	Often it is useful to look at the frequency response on
%	logarithmic (dB) scales.

pause % Strike any key for plot.

plot(ff(2:n),20*log10(h(2:n,:))), ...
title('YuleWalk, Butterworth and Chebyshev filters'), ...
xlabel('Frequency (Hz)'), ylabel('Magnitude in dB'),  pause
clc
