|
wolffd@0
|
1 function y=rasta(x,fs,low,high)
|
|
wolffd@0
|
2 % function y=rasta(x,fs) where x is the input data (rows of time data),
|
|
wolffd@0
|
3 % and fs is the frame rate (sampling rate) in Hz. This is a modified
|
|
wolffd@0
|
4 % version of the original filter. Here the RASTA filter is approximated
|
|
wolffd@0
|
5 % by a simple fourth order Butterworth bandpass filter. See pages 50-51
|
|
wolffd@0
|
6 % of my second IRC logbook for the derivation.
|
|
wolffd@0
|
7 %
|
|
wolffd@0
|
8 % Hermansky and Morgan, "RASTA Processing of Speech." IEEE Transactions
|
|
wolffd@0
|
9 % on Speech and Audio Processing. vol. 2, no. 4, October 1994
|
|
wolffd@0
|
10 %
|
|
wolffd@0
|
11
|
|
wolffd@0
|
12 % (c) 1998 Interval Research Corporation
|
|
wolffd@0
|
13 % Malcolm Slaney, January 30, 1996, Interval Research Corporation
|
|
wolffd@0
|
14
|
|
wolffd@0
|
15 if (nargin < 2); fs=100; end
|
|
wolffd@0
|
16 if (nargin < 3); low=.9; end
|
|
wolffd@0
|
17 if (nargin < 4); high=12.8; end
|
|
wolffd@0
|
18
|
|
wolffd@0
|
19 if (low == 0 & high == 0) % Original Filter
|
|
wolffd@0
|
20 num = .1*[2 1 0 -1 -2];
|
|
wolffd@0
|
21 denum = [1 -.94];
|
|
wolffd@0
|
22 else % New fourth order
|
|
wolffd@0
|
23 % Butterworth BP filter
|
|
wolffd@0
|
24 w1=low/fs*2*pi;
|
|
wolffd@0
|
25 w2=high/fs*2*pi;
|
|
wolffd@0
|
26 theta=1;
|
|
wolffd@0
|
27
|
|
wolffd@0
|
28 a=cos((w1+w2)/2)/cos((w2-w1)/2);
|
|
wolffd@0
|
29 k=cot((w2-w1)/2)*tan(theta/2);
|
|
wolffd@0
|
30
|
|
wolffd@0
|
31 num = [1 0 -2 0 1];
|
|
wolffd@0
|
32 denum = [(1 + 2*2^(1/2)*k + 4*k^2) ...
|
|
wolffd@0
|
33 (-4*2^(1/2)*a*k - 16*a*k^2) ...
|
|
wolffd@0
|
34 (-2 + 8*k^2 + 16*a^2*k^2) ...
|
|
wolffd@0
|
35 (4*2^(1/2)*a*k - 16*a*k^2) ...
|
|
wolffd@0
|
36 (1 - 2*2^(1/2)*k + 4*k^2)];
|
|
wolffd@0
|
37 scale = denum(1); % Scale by a(1) component
|
|
wolffd@0
|
38 num = num/scale;
|
|
wolffd@0
|
39 denum = denum/scale;
|
|
wolffd@0
|
40 end
|
|
wolffd@0
|
41
|
|
wolffd@0
|
42 if (0)
|
|
wolffd@0
|
43 len = 1024;
|
|
wolffd@0
|
44 impulse = zeros(1,len);
|
|
wolffd@0
|
45 impulse(1) = 1;
|
|
wolffd@0
|
46
|
|
wolffd@0
|
47 y=filter(num,denum,impulse);
|
|
wolffd@0
|
48 ym = abs(fft(y));
|
|
wolffd@0
|
49 ym=20*log10(ym);
|
|
wolffd@0
|
50 f=(0:(len-1))/len*fs;
|
|
wolffd@0
|
51 semilogx(f(1:len/2),ym(1:len/2));
|
|
wolffd@0
|
52 drawnow;
|
|
wolffd@0
|
53 end
|
|
wolffd@0
|
54
|
|
wolffd@0
|
55 if (length(x) == size(x,1)*size(x,2))
|
|
wolffd@0
|
56 y = filter(num,denum,x,x(1)*[-1 -1 1 1]);
|
|
wolffd@0
|
57 else
|
|
wolffd@0
|
58 y = zeros(size(x));
|
|
wolffd@0
|
59 for i=1:size(x,1)
|
|
wolffd@0
|
60 y(i,:) = filter(num,denum,x(i,:),x(i,1)*[-1 -1 1 1]/scale);
|
|
wolffd@0
|
61 end
|
|
wolffd@0
|
62 end
|
