fc = 1 ; type = 'LP' ; fs = 10*fc ;
%
%	Lowpass, Rectangular Window
%
h = FIRDF(type,fc,fs,'RE',0.01) ;	% Return impulse response of filter
N = length(h) ; Ts = 1/fs ; n = [0:N-1]' ;
F = [0:0.001:1/2]' ; 	% Vector of DT frequencies for plotting frequency response

[H,F] = DTFT(n,h,F) ;	% Find the DTFT of the impulse response
subplot(2,2,1) ;
p = xyplot(F*fs,20*log10(abs(H)),[0,fs/2,-120,0],'\itf ',...
			'|H(\ite^{{\itj}2{\pi}{\itf}{\itT}_s} )|',...
			'Times',18,'Times',14,...
'Lowpass - Rectangular Window','Times',24,'n','c') ;
subplot(2,2,2) ;
p = xyplot(F*fs,20*log10(abs(H)),[0,fs/10,-5,0],'\itf ',...
			'|H(\ite^{{\itj}2{\pi}{\itf}{\itT}_s} )|',...
			'Times',18,'Times',14,...
'Lowpass - Rectangular Window','Times',24,'n','c') ;
%
%	Lowpass, von Hann Window
%
h = FIRDF(type,fc,fs,'VH',0.01) ;
N = length(h) ; Ts = 1/fs ; n = [0:N-1]' ;
F = [0:0.001:1/2]' ; [H,F] = DTFT(n,h,F) ;
subplot(2,2,3) ;
p = xyplot(F*fs,20*log10(abs(H)),[0,fs/2,-120,20],'\itf ',...
			'|H(\ite^{{\itj}2{\pi}{\itf}{\itT}_s} )|',...
			'Times',18,'Times',14,...
'Lowpass - Von Hann Window','Times',24,'n','c') ;
subplot(2,2,4) ;
p = xyplot(F*fs,20*log10(abs(H)),[0,fs/10,-5,0],'\itf ',...
			'|H(\ite^{{\itj}2{\pi}{\itf}{\itT}_s} )|',...
			'Times',18,'Times',14,...
'Lowpass - Von Hann Window','Times',24,'n','c') ;


%	Function to design a digital filter using truncation of the 
%	impulse response to approximate the ideal impulse response.  The 
%	user specifies the type of ideal filter,
%
%	LP	-	lowpass
%	BP	-	bandpass
%
%	the cutoff frequency(s) fcs (a scalar for LP and HP a 2-vector
%	for BP and BS), the sampling rate, fs, the type of window,
%
%	RE	-	rectangular
%	VH	-	von Hann
%	BA	-	Bartlett
%	HA	-	Hamming
%	BL	-	Blackman
%
%	and the allowable truncation error,
%	err, as a fraction of the total impulse response signal energy.
%
%	The function returns the filter coefficients as a vector.
%
%	function h = FIRDF(type,fcs,fs,window,err)
%
function h = FIRDF(type,fcs,fs,window,err)

	if fs == 0 | fs == inf | fs == -inf,
		disp('Sampling rate is unusable') ;
	else
		fs
		Ts = 1/fs ; 
		type = upper(type) ; window = upper(window) ;
		switch type
			case 'LP',	%  Lowpass filter design
				fc = fcs(1) ; 
				zc1 = 1/(2*fc) ;	% "zc" stands for zero-crossing
				N = 200*zc1/Ts ; n = [0:N]' ;
				Etotal = sum(sinc(2*fc*n*Ts).^2) ;	% Total signal energy
				E = 0 ; N = 1 ;
				while abs(E-Etotal) > Etotal*err,	% Find duration required for
									% error specified
					N = 2*N ;
					n = [0:N]' ;
					E = sum(sinc(2*fc*n*Ts).^2) ;
				end
				delN = floor(N/10) ;
				while abs(E-Etotal) < Etotal*err & abs(delN) > 1,
					N = N-delN ;
					n = [0:N]' ;
					E = sum(sinc(2*fc*n*Ts).^2) ;
					delN = floor(delN/2) ;
				end
				N = ceil(N/2)*2 ;	%	Make number of pts even
				nmid = (N/2-1) ;
				w = makeWindow(window,n,N) ;
				h = w.*sinc(2*fc*(n-nmid)*Ts) ;	% Truncated impulse response
				h = h/sum(h) ;		% Set zero-frequency gain to one

			case 'BP'	% Bandpass filter design
				fl = fcs(1) ; fh = fcs(2) ; 	% Low and high cutoff frequencies
				fmid = (fh + fl)/2 ; 		% Center frequency
				df = abs(fh-fl) ; 		% Filter bandwidth
				zc1 = 1/df ;			% "zc" stands for zero crossing
				N = 200*zc1/Ts ; n = [0:N]' ;
				Etotal = sum((2*df*sinc(df*n*Ts).*...
						cos(2*pi*fmid*n*Ts)).^2) ;	% Total signal energy
				E = 0 ; N = 1 ;
				while abs(E-Etotal) > Etotal*err,	% Find the duration for the specified
									% error
					N = 2*N ;
					n = [0:N]' ;
					E = sum((2*df*sinc(df*n*Ts).*...
						cos(2*pi*fmid*n*Ts)).^2) ;
				end
				delN = floor(N/10) ;
				while abs(E-Etotal)<Etotal*err & abs(delN)>1,
					N = N-delN ;
					n = [0:N]' ;
					E = sum((2*df*sinc(df*n*Ts).*...
						cos(2*pi*fmid*n*Ts)).^2) ;
					delN = floor(delN/2) ;
				end
				N = ceil(N/2)*2 ;	%	Make number of pts even
				nmid = (N/2-1) ; n = [0:N-1]' ;
				w = makeWindow(window,n,N) ;
				h = w.*2.*df.*sinc(df*(n-nmid)*Ts).*...
						cos(2*pi*fmid*(n-nmid)*Ts) ;	% Truncated impulse response
				h = h/sum(h.*cos(2*pi*fmid*(n-nmid)*Ts)) ;	% Set zero-frequency gain to one
		end
	end

function w = makeWindow(window,n,N)
	switch window
		case 'RE'
			w = ones(N,1) ;
		case 'VH'
			w = (1-cos(2*pi*n/(N-1)))/2 ;
		case 'BA'
			w = 2*n/(N-1).*(0<=n & n<=(N-1)/2) + ...
					(2-2*n/(N-1)).*((N-1)/2<=n & n<N) ;
		case 'HA'
			w = 0.54 - 0.46*cos(2*pi*n/(N-1)) ;
		case 'BL'
			w = 0.42 - 0.5*cos(2*pi*n/(N-1)) + ...
						0.08*cos(4*pi*n/(N-1)) ;
		otherwise
			w = ones(N,1) ;
	end

%	Function to compute an approximation to the discrete-time Fourier transform, 
%	X(F), of a set of input data, x[n], n = 0 to NF-1 at the discrete-time
%	frequencies, F.  If the input vector of DT frequencies, F, is omitted, the 
%	function returns the values of X(F) at the harmonics of the fundamental
%	discrete-time frequency, 0, 1/NF, 2/NF,...(NF-1)/NF.
%
%	n and x must be column vectors of real numbers and must have the same length
%	F must be a column vector of real numbers
%
%	[X,F] = DTFT(n,x,F)

function  [X,F] = DTFT(n,x,F)
	if nargin == 3,
		Fvec = [0:length(x)-1]'/length(x) ;
		Xper = DFT(x) ;
		Xper = [Xper ; Xper(1)] ; Fvec = [Fvec ; 1] ;
% 		Xper = Xper.*exp(-j*2*pi*n(1)*Fvec) ;
		X = interp1(Fvec,Xper,mod(F,1),'cubic') ;
	else
		F = [0:length(x)-1]'/length(x) ; X = DFT(x).*exp(-j*2*pi*n(1)*F) ;
	end
	
%	Function to compute the discrete Fourier transform, X[k], of an input set of 
%	data, x[n], n = 0 to NF-1, at the harmonic numbers, k.  The vector returned 
%	is a vector of X[k]'s at the k's input to the function.  If only the "x" 
%	argument is sent, the function returns the NF data for the harmonic-number 
%	period, k = 0 to NF-1.
%
%	x must be a column vector of real numbers
%
%	[X,k] = DFT(x,k)

function [X,k] = DFT(x,k)
	if nargin == 2,
		Xper = fft(x) ; X = Xper(mod(k,length(x))+1) ;
	else
		k = (0:length(x)-1)' ; X = fft(x) ;
	end

