comparison toolboxes/MIRtoolbox1.3.2/AuditoryToolbox/rasta.m @ 0:e9a9cd732c1e tip

first hg version after svn

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