filtbank.m

来自「车辆工程平顺性方面的1/3倍中心频率的有关程序」· M 代码 · 共 76 行

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));
end


if (Fc(i_low) > Fs/20)
  i_dec = 0; 
else
  i_dec = max(find(Fc<=Fs/20)); 
end


for i = i_up:-1:i_dec+1
  [B,A] = oct3dsgn(Fc(i),Fs,N);
  y = filter(B,A,x); 
  P(:,i) = leq(y,T);
end


if (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);    
  end
end

f = Fref(i_low:i_up);
p = P(:,i_low:i_up);