wolffd@0: function y=rasta(x,fs,low,high)
wolffd@0: % function y=rasta(x,fs) where x is the input data (rows of time data), 
wolffd@0: % and fs is the frame rate (sampling rate) in Hz.  This is a modified 
wolffd@0: % version of the original filter. Here the RASTA filter is approximated
wolffd@0: % by a simple fourth order Butterworth bandpass filter.  See pages 50-51
wolffd@0: % of my second IRC logbook for the derivation.
wolffd@0: %
wolffd@0: % Hermansky and Morgan, "RASTA Processing of Speech."  IEEE Transactions
wolffd@0: % on Speech and Audio Processing.  vol. 2, no. 4, October 1994
wolffd@0: %
wolffd@0: 
wolffd@0: % (c) 1998 Interval Research Corporation  
wolffd@0: % Malcolm Slaney, January 30, 1996, Interval Research Corporation
wolffd@0: 
wolffd@0: if (nargin < 2); fs=100; end
wolffd@0: if (nargin < 3); low=.9; end
wolffd@0: if (nargin < 4); high=12.8; end
wolffd@0: 
wolffd@0: if (low == 0 & high == 0)			% Original Filter
wolffd@0: 	num = .1*[2 1 0 -1 -2];
wolffd@0: 	denum = [1 -.94];
wolffd@0: else						% New fourth order 
wolffd@0: 						% Butterworth BP filter
wolffd@0: 	w1=low/fs*2*pi;
wolffd@0: 	w2=high/fs*2*pi;
wolffd@0: 	theta=1;
wolffd@0: 	
wolffd@0: 	a=cos((w1+w2)/2)/cos((w2-w1)/2);
wolffd@0: 	k=cot((w2-w1)/2)*tan(theta/2);
wolffd@0: 	
wolffd@0: 	num = [1 0 -2 0 1];
wolffd@0: 	denum = [(1 + 2*2^(1/2)*k + 4*k^2) ...
wolffd@0: 			(-4*2^(1/2)*a*k - 16*a*k^2) ...
wolffd@0: 			(-2 + 8*k^2 + 16*a^2*k^2) ...
wolffd@0: 			(4*2^(1/2)*a*k - 16*a*k^2) ...
wolffd@0: 			(1 - 2*2^(1/2)*k + 4*k^2)];
wolffd@0: 	scale = denum(1);			% Scale by a(1) component
wolffd@0: 	num = num/scale;
wolffd@0: 	denum = denum/scale;
wolffd@0: end
wolffd@0: 
wolffd@0: if (0)
wolffd@0: 	len = 1024;
wolffd@0: 	impulse = zeros(1,len);
wolffd@0: 	impulse(1) = 1;
wolffd@0: 	
wolffd@0: 	y=filter(num,denum,impulse);
wolffd@0: 	ym = abs(fft(y));
wolffd@0: 	ym=20*log10(ym);
wolffd@0: 	f=(0:(len-1))/len*fs;
wolffd@0: 	semilogx(f(1:len/2),ym(1:len/2));
wolffd@0: 	drawnow;
wolffd@0: end
wolffd@0: 
wolffd@0: if (length(x) == size(x,1)*size(x,2))
wolffd@0: 	y = filter(num,denum,x,x(1)*[-1 -1 1 1]);
wolffd@0: else
wolffd@0: 	y = zeros(size(x));
wolffd@0: 	for i=1:size(x,1)
wolffd@0: 		y(i,:) = filter(num,denum,x(i,:),x(i,1)*[-1 -1 1 1]/scale);
wolffd@0: 	end
wolffd@0: end