tomwalters@0: %
tomwalters@0: %	Gammachirp : Theoretical auditory filter 
tomwalters@0: %	Toshio IRINO
tomwalters@0: %	7 Apr. 97 (additional comments)
tomwalters@0: %	20 Aug. 97 (Simplify & Carrier Selection)
tomwalters@0: %	10 Jun. 98 (SwNorm)
tomwalters@0: %	26 Nov. 98 (phase = phase + c ln fr/f0)
tomwalters@0: %	7  Jan. 2002 (adding 'envelope' option)
tomwalters@0: %	22  Nov. 2002 (debugging 'peak' option)
tomwalters@0: %
tomwalters@0: %	gc(t) = t^(n-1) exp(-2 pi b ERB(Frs)) cos(2*pi*Frs*t + c ln t + phase)
tomwalters@0: %
tomwalters@0: %	function [GC, LenGC, Fps, InstFreq ] ...
tomwalters@0: %	       = GammaChirp(Frs,SR,OrderG,CoefERBw,CoefC,Phase,SwCarr,SwNorm);
tomwalters@0: %	INPUT : Frs	: Asymptotic Frequency ( vector )
tomwalters@0: %		SR 	: Sampling Frequency
tomwalters@0: %		OrderG 	: Order of Gamma function t^(OrderG-1)        == n   
tomwalters@0: %		CoefERBw: Coeficient -> exp(-2*pi*CoefERBw*ERB(f))    == b
tomwalters@0: %		CoefC	: Coeficient -> exp(j*2*pi*Frs + CoefC*ln(t)) == c
tomwalters@0: %		Phase	: Start Phase(0 ~ 2*pi)                       
tomwalters@0: %		SwCarr  : Carrier ('cos','sin','complex','envelope': 3 letters)
tomwalters@0: %		SwNorm  : Normalization of peak spectrum level ('no', 'peak')
tomwalters@0: %	OUTPUT: GC 	: GammaChirp                     ( matrix )
tomwalters@0: %		LenGC 	: Length of GC for each channel  ( vector )
tomwalters@0: %               Fps     : Peak Frequency                 ( vector )
tomwalters@0: %		InstFreq: Instanteneous Frequency        ( matrix )
tomwalters@0: %
tomwalters@0: %	
tomwalters@0: function [GC, LenGC, Fps, InstFreq ] ...
tomwalters@0:     = GammaChirp(Frs,SR,OrderG,CoefERBw,CoefC,Phase,SwCarr,SwNorm);
tomwalters@0: 
tomwalters@0: if nargin < 2,            help GammaChirp; return; end;
tomwalters@0: Frs = Frs(:);
tomwalters@0: NumCh = length(Frs);
tomwalters@0: if nargin < 3,            OrderG = [];	               end;
tomwalters@0: if length(OrderG) == 0,   OrderG = 4;                  end; % Default GammaTone
tomwalters@0: if length(OrderG) == 1,   OrderG = OrderG*ones(NumCh,1); end;
tomwalters@0: if nargin < 4,            CoefERBw = [];	       end; 
tomwalters@0: if length(CoefERBw) == 0, CoefERBw = 1.019;            end; % Default GammaTone
tomwalters@0: if length(CoefERBw) == 1, CoefERBw = CoefERBw*ones(NumCh,1); end; 
tomwalters@0: if nargin < 5,            CoefC  = [];                 end;
tomwalters@0: if length(CoefC) == 0,    CoefC  = 0;                  end; % Default GammaTone
tomwalters@0: if length(CoefC) == 1,    CoefC  = CoefC*ones(NumCh,1); end; 
tomwalters@0: if nargin < 6,            Phase  = [];                  end;
tomwalters@0: if length(Phase) == 0,    Phase  = 0;                  end; 
tomwalters@0: if length(Phase) == 1,    Phase  = Phase*ones(NumCh,1); end; 
tomwalters@0: if nargin < 7,            SwCarr = [];                 end;
tomwalters@0: if length(SwCarr) == 0,   SwCarr = 'cos';	       end;
tomwalters@0: if nargin < 8,            SwNorm = [];                 end;
tomwalters@0: if length(SwNorm) == 0,   SwNorm = 'no';	       end;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: [ERBrate ERBw] = Freq2ERB(Frs);                             % G&M (1990)
tomwalters@0: LenGC1kHz = (40*max(OrderG)/max(CoefERBw) + 200)*SR/16000;  % 2 Aug 96 
tomwalters@0: [dummy ERBw1kHz] = Freq2ERB(1000);	
tomwalters@0: 
tomwalters@0: if strcmp(SwCarr,'sin'), Phase = Phase - pi/2*ones(1,NumCh); end;
tomwalters@0: %%% Phase compensation
tomwalters@0: Phase = Phase + CoefC.*log(Frs/1000); % relative phase to 1kHz
tomwalters@0: 
tomwalters@0: LenGC = fix(LenGC1kHz*ERBw1kHz./ERBw);
tomwalters@0: 
tomwalters@0: %%%%%  Production of GammaChirp  %%%%%
tomwalters@0: GC       = zeros(NumCh,max(LenGC));
tomwalters@0: if nargout > 2, Fps = Fr2Fpeak(OrderG,CoefERBw,CoefC,Frs); end; % Peak Freq.
tomwalters@0: if nargout > 3, InstFreq = zeros(NumCh,max(LenGC));        end;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: for nch = 1:NumCh,
tomwalters@0:   t = (1:LenGC(nch)-1)/SR;
tomwalters@0: 
tomwalters@0:   GammaEnv = t.^(OrderG(nch)-1).*exp(-2*pi*CoefERBw(nch)*ERBw(nch)*t);
tomwalters@0:   GammaEnv = [ 0 GammaEnv/max(GammaEnv)];
tomwalters@0: 
tomwalters@0:   if strcmp(SwCarr(1:3),'env') % envelope
tomwalters@0:     Carrier = ones(size(GammaEnv));
tomwalters@0:   elseif strcmp(SwCarr(1:3),'com') % complex
tomwalters@0:     Carrier = [ 0 exp(i * (2*pi*Frs(nch)*t + CoefC(nch)*log(t) +Phase(nch)) )];
tomwalters@0:   else
tomwalters@0:     Carrier = [ 0 cos(2*pi*Frs(nch)*t + CoefC(nch)*log(t) +Phase(nch))];
tomwalters@0:   end;
tomwalters@0: 
tomwalters@0:   GC(nch,1:LenGC(nch)) = GammaEnv.*Carrier;
tomwalters@0: 
tomwalters@0:   if nargout > 3, 
tomwalters@0:     InstFreq(nch,1:LenGC(nch)) = [0, [Frs(nch) + CoefC(nch)./(2*pi*t)]];
tomwalters@0:   end;
tomwalters@0: 
tomwalters@0:   if strcmp(SwNorm,'peak') == 1,  % peak gain normalization
tomwalters@0:      [frsp freq] = freqz(GC(nch,1:LenGC(nch)),1,LenGC(nch),SR);
tomwalters@0:      fp = Fr2Fpeak(OrderG(nch),CoefERBw(nch),CoefC(nch),Frs(nch));
tomwalters@0:      [dummy np] = min(abs(freq-fp)); 
tomwalters@0:      GC(nch,:) = GC(nch,:)/abs(frsp(np));
tomwalters@0:   end;
tomwalters@0: 
tomwalters@0: end; % nch = ...
tomwalters@0: 
tomwalters@0: return
tomwalters@0:     
tomwalters@0: %% ERBw = 0.128*Frs;     % Complete Constant Q only for check.
tomwalters@0: 
tomwalters@0: % old 
tomwalters@0: % Amp = ones(NumCh,1);                                  % No normalization
tomwalters@0: % if strcmp(SwNorm,'peak'),  Amp = ERBw./ERBw1kHz; end; % Peak spectrum==const. 
tomwalters@0: % when it is gammatone
tomwalters@0: %  if strcmp(SwNorm,'peak'),   ...
tomwalters@0: %         Amp = 2.815*sqrt(4/OrderG).*CoefERBw.*ERBw/SR; end;
tomwalters@0: % Peak spectrum==const. The gain is 1.0 when filtering sinusoid at cf.
tomwalters@0: % GC(nch,:) = GC(nch,:)/max(abs(freqz(GC(nch,:),1,LenGC(nch))));
tomwalters@0: %