% COMBFFT	Graphically shows the relationship of the FFT to processing a
%		time window of data through a sequence of "COMB Filter" banks.
%		It will sequentially display each filter bank's frequency response
%		then show the effective	concatenated filter. 
%		This process is repeated for each of the 16 output "bins" of the
%		16 point, decimated-in-time, N=16 FFT.

clear
clg
clc
Y = zeros(1,128);

for k = 0:1:15

	%  Defines the phase factor associated with
	%  each filter bank in the form of 
	%  y(k) = x(k) + exp(j*theta*N)*x(k-N)
	%  to produce each term of the FFT.

	disp('Press the spacebar to continue when it pauses')

	T1 = exp(j*k*pi/1);	
	T2 = exp(j*k*pi/2);
	T3 = exp(j*k*pi/4);
	T4 = exp(j*k*pi/8);	
	
	% H(z) for first stage of FFT
	num = [1 0 0 0 0 0 0 0 T1];
	den = [2 0 0 0 0 0 0 0 0 ];	% the '2' normalizes the output to 
							% unity gain for this stage
	[H1,q] = freqz(num,den,128,'whole');
			
	% H(z) for second stage of FFT
	num = [1 0 0 0 T2];
	den = [2 0 0 0 0];			% the '2' normalizes the output to 
							% unity gain for this stage
	[H2,q] = freqz(num,den,128,'whole');

	% H(z) for third stage of FFT
	num = [1 0 T3];
	den = [2 0 0];				% the '2' normalizes the output to 
							% unity gain for this stage
	[H3,q] = freqz(num,den,128,'whole');

	% H(z) for fourth stage of FFT	
	num = [1 T4];
	den = [2 0];				% the '2' normalizes the output to 
							% unity gain for this stage
	[H4,q] = freqz(num,den,128,'whole');
		
	% H(z) for output of FFT
	H5 = H1 .* H2 .* H3 .* H4;

	% Frequency resolution term of DFT/FFT
	q = (q*8)/pi;
	
	figure(1)
	plot(q,abs(H1))
	grid
	xlabel('Frequency resolution in terms of N')
	ylabel('Magnitude')
	title('First Stage')
	pause
	
	figure(2)
	plot(q,abs(H2))
	grid
	xlabel('Frequency resolution in terms of N')
	ylabel('Magnitude')
	title('Second Stage')
     pause	
	
	figure(3)
	plot(q,abs(H3))
	grid
	xlabel('Frequency resolution in terms of N')
	ylabel('Magnitude')
	title('Third Stage')
     pause	
	
	figure(4)
	plot(q,abs(H4))
	grid
	xlabel('Frequency resolution in terms of N')
	ylabel('Magnitude')
	title('Fourth Stage')	
     pause	
	
	figure(5)
	plot(q,abs(H5))
	grid
	xlabel('Frequency resolution in terms of N')
	ylabel('Magnitude')
	title(['Output for k = ' num2str(k)])
	pause	
	
	% Test to see that [a] passed through comb filter
	% banks will result in FFT/DFT.  Uses convolution
	% property of fourier transforms in this stagement.
	% P(omega) = A(omega)*H(omega) where H is the
	% frequency response of the comb filter, and
	% A is the frequency response of the 8 point input
	% sequence.
	% a = [1 0 0 0.5 0 0 1 0 0 0 0 0 1 0 0 0];
	% A = fft(a,128);
	% P = A.' .* H5;
	% p = ifft(P);

	% Y(k+1) = p(1);
	
end

% Table listing actual fft/dft of [a] and Y is the resulting
% output of comb filter sequence.
% [fft(a); Y]