tomwalters@0: %
tomwalters@0: %	Frequency Response of GammaChirp
tomwalters@0: %	Toshio IRINO
tomwalters@0: %	Original Version: 30 Sept 98
tomwalters@0: %	Modified:          7 June 2004 (removing transpose when NumCh==1)
tomwalters@0: %
tomwalters@0: %    function [AmpFrsp, freq, Fpeak, GrpDly, PhsFrsp] ...
tomwalters@0: %       = GammaChirpFrsp(Frs,SR,OrderG,CoefERBw,CoefC,Phase,NfrqRsl);
tomwalters@0: %	INPUT : Frs	: Resonance Freq. (vector)
tomwalters@0: %		SR 	: Sampling Freq.
tomwalters@0: %		OrderG 	: Order of Gamma function t^(OrderG-1) (vector)
tomwalters@0: %		CoefERBw: Coeficient -> exp(-2*pi*CoefERBw*ERB(f))
tomwalters@0: %		CoefC	: Coeficient -> exp(j*2*pi*Fr + CoefC*ln(t))
tomwalters@0: %		Phase	: Start Phase (0 ~ 2*pi)
tomwalters@0: %		NfrqRsl : freq. resolution 
tomwalters@0: %	OUTPUT: AmpFrsp	: abs(Response)      (NumCh*NfrqRsl matrix)
tomwalters@0: %		freq	: frequency          (1 * NfrqRsl vector)
tomwalters@0: %		Fpeak	: Peak frequency     (NumCh * 1 vector)
tomwalters@0: %               GrpDly  : Group delay        (NumCh*NfrqRsl matrix)
tomwalters@0: %               PhsFrsp : angle(Response);   (NumCh*NfrqRsl matrix)
tomwalters@0: %
tomwalters@0: function [AmpFrsp, freq, Fpeak, GrpDly, PhsFrsp ] ...
tomwalters@0:     = GammaChirpFrsp(Frs,SR,OrderG,CoefERBw,CoefC,Phase,NfrqRsl);
tomwalters@0: 
tomwalters@0: if nargin < 1, help GammaChirpFrsp; return; end;
tomwalters@0: if nargin < 2, SR = 48000; 	end;
tomwalters@0: if length(SR) == 0, error('Specify Sampling Frequency'); end;
tomwalters@0: Frs = Frs(:);
tomwalters@0: NumCh = length(Frs);
tomwalters@0: 
tomwalters@0: if nargin < 3, OrderG	= 4;                    end; 	% Default GammaTone
tomwalters@0: if length(OrderG)   == 1, OrderG = OrderG*ones(NumCh,1); end;
tomwalters@0: if nargin < 4, CoefERBw = 1.019;	        end; 	% Default GammaTone
tomwalters@0: if length(CoefERBw) == 1, CoefERBw = CoefERBw*ones(NumCh,1); end;
tomwalters@0: if nargin < 5, 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 = zeros(NumCh,1);	end;	% Default GammaTone
tomwalters@0: if nargin < 7, NfrqRsl  = 1024;                end;
tomwalters@0: if NfrqRsl < 256, help GammaChirpFrsp;  error('NfrqRsl < 256'); end;
tomwalters@0: 
tomwalters@0: [ERBrate ERBw] = Freq2ERB(Frs(:));
tomwalters@0: freq = [0:NfrqRsl-1]/NfrqRsl*SR/2;
tomwalters@0: 
tomwalters@0: one1 = ones(1,NfrqRsl);
tomwalters@0: bh = ( CoefERBw(:).*ERBw(:) ) * one1;
tomwalters@0: fd = ones(NumCh,1)*freq(:)' - Frs(:)*one1;
tomwalters@0: cn = ( CoefC(:)./OrderG(:) ) * one1;
tomwalters@0: n  = OrderG(:) *one1;
tomwalters@0: c  = CoefC(:)  *one1;
tomwalters@0: Phase = Phase(:) *one1;
tomwalters@0: 
tomwalters@0: %%% Analytic form (normalized at Fpeak) %%%
tomwalters@0: AmpFrsp = ( (1 + cn.^2)./(1 + (fd./bh).^2) ).^(n/2) ...
tomwalters@0:       .* exp(  c .* (atan(fd./bh) - atan(cn)));
tomwalters@0: 
tomwalters@0: if nargout > 2,
tomwalters@0:   Fpeak = Frs + CoefERBw.*ERBw.*CoefC./OrderG;
tomwalters@0:   if nargout > 3,
tomwalters@0:     GrpDly = 1/(2*pi)*(n.*bh + c.*fd)./(bh.^2 + fd.^2);
tomwalters@0:     if nargout > 4,
tomwalters@0:       PhsFrsp= -n.*atan(fd./bh) -c./2.*log((2*pi*bh).^2 +(2*pi*fd).^2 )+Phase;
tomwalters@0:     end;
tomwalters@0:   end;
tomwalters@0: end;
tomwalters@0: 
tomwalters@0: return
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: %%% No use anymore since it is somewhat confusing  7 June 2004 %%%%%%%%%%%%
tomwalters@0: 
tomwalters@0: if NumCh == 1, % to let the results in column vector (same as freqz)
tomwalters@0:   AmpFrsp = AmpFrsp(:); 
tomwalters@0:   freq=freq(:); 
tomwalters@0:   if nargout > 3, 
tomwalters@0:     GrpDly =  GrpDly(:);
tomwalters@0:     if nargout > 4,
tomwalters@0:       PhsFrsp =  PhsFrsp(:);
tomwalters@0:     end;
tomwalters@0:   end;
tomwalters@0: end; 
tomwalters@0: 
tomwalters@0: return
tomwalters@0: 
tomwalters@0: 
tomwalters@0: % Fbw   = sqrt( -1 + 2.^(1./OrderG)).*CoefERBw.*ERBw; % bandwidth when c=0
tomwalters@0: 
tomwalters@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tomwalters@0: 
tomwalters@0: % OLD
tomwalters@0: % val = 2*pi*(bh + i*fd);
tomwalters@0: % arg = angle(val);
tomwalters@0: % FrspAna = 1./(abs(val)).^OrderG(nch) .* exp(CoefC(nch)* arg);
tomwalters@0: % FrspAna = FrspAna/max(FrspAna);
tomwalters@0:   
tomwalters@0: 
tomwalters@0: if 0
tomwalters@0:   
tomwalters@0: for nch = 1:NumCh
tomwalters@0:   bh = CoefERBw(nch)*ERBw(nch);
tomwalters@0:   fd = freq(:)' - Frs(nch);
tomwalters@0:   cn = CoefC(nch)/OrderG(nch);
tomwalters@0:   %%% Analytic form (normalized at Fpeak) %%%
tomwalters@0:   FrspAna = (sqrt( (1+cn^2)./(1+(fd/bh).^2) ) ).^OrderG(nch) ...
tomwalters@0:       .* exp(  CoefC(nch) .* (atan(fd/bh) - atan(cn)));
tomwalters@0:   %  Frsp(nch, 1:NfrqRsl) = FrspAna;
tomwalters@0:   sqrt(mean( ( Frsp(nch, 1:NfrqRsl) - FrspAna ).^2) )
tomwalters@0:   GrpDly(nch,1:NfrqRsl) = ...
tomwalters@0:       1/(2*pi)*(OrderG(nch)*bh + CoefC(nch)*fd)./(bh.^2 + fd.^2);
tomwalters@0: end;
tomwalters@0: 
tomwalters@0: end;
tomwalters@0: 
tomwalters@0: