function [p,f] = filtbank(x,Fs)Fref = [1.6 2 2.5,3.15 4 5,6.3 8 10,12.5 16 20,25 31.5 40,50 63 80];        Fc = 1.6*((2^(1/3)).^[0:1:17]);       N = 3;  				  U = 2^(1/3); T = length(x); P = zeros(1,length(Fref));i_up = 18; 				i_low = 1;if (Fs/2) < Fref(i_low)*1.5   error('Sampling frequency Fs too low.'); elseif (Fs/2) < Fref(i_up)*1.5  disp('Warning: frequency range must be reduced (Fs too low).');   i_up = max(find(Fc<=Fs/3));endif (Fc(i_low) > Fs/20)  i_dec = 0; else  i_dec = max(find(Fc<=Fs/20)); endfor i = i_up:-1:i_dec+1  [B,A] = oct3dsgn(Fc(i),Fs,N);  y = filter(B,A,x);   P(:,i) = leq(y,T);endif (i_dec > 0)  [Bu,Au] = oct3dsgn(Fc(i_dec),Fs/2,N);   [Bc,Ac] = oct3dsgn(Fc(i_dec)/U,Fs/2,N);   [Bl,Al] = oct3dsgn(Fc(i_dec)/(U^2),Fs/2,N);   i = i_dec;   while  i >= i_low+2     x = decimate(x,2);      T = T/2;     y = filter(Bu,Au,x);      P(:,i) =  leq(y,T);         y = filter(Bc,Ac,x);      P(:,i-1) =  leq(y,T);         y = filter(Bl,Al,x);      P(:,i-2) =  leq(y,T);      i = i-3;   end   if (i == (i_low+1))     x = decimate(x,2);      T = T/2;     y = filter(Bu,Au,x);      P(:,i) =  leq(y,T);         y = filter(Bc,Ac,x);      P(:,i-1) =  leq(y,T);      elseif (i == (i_low))     x = decimate(x,2);      T = T/2;     y = filter(Bu,Au,x);      P(:,i) =  leq(y,T);      endendf = Fref(i_low:i_up);p = P(:,i_low:i_up);