[5080e1]: auditory / gammatonefir.m Maximize Restore History

Download this file

gammatonefir.m    135 lines (107 with data), 4.1 kB

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
function b=gammatonefir(fc,fs,varargin);
%GAMMATONEFIR Gammatone filter coefficients
% Usage: b = gammatonefir(fc,fs,n,betamul);
% b = gammatonefir(fc,fs,n);
% b = gammatonefir(fc,fs);
%
% Input parameters:
% fc : center frequency in Hz.
% fs : sampling rate in Hz.
% n : filter order.
% beta : bandwidth of the filter.
%
% Output parameters:
% b : FIR filters as columns
%
% `gammatonefir(fc,fs,n,betamul)` computes the filter coefficients of a
% digital FIR gammatone filter of length *n* with center frequency *fc*,
% 4th order rising slope, sampling rate *fs* and bandwith determined by
% *betamul*. The bandwidth *beta* of each filter is determined as *betamul*
% times |audfiltbw|_ of the center frequency of corresponding filter.
%
% `gammatonefir(fc,fs,n)` will do the same but choose a filter bandwidth
% according to Glasberg and Moore (1990). *betamul* is choosen to be 1.0183.
%
% `gammatonefir(fc,fs)` will do as above and choose a sufficiently long
% filter to accurately represent the lowest subband channel.
%
% If *fc* is a vector, each entry of *fc* is considered as one center
% frequency, and the corresponding coefficients are returned as column
% vectors in the output.
%
% The inpulse response of the gammatone filter is given by
%
% .. g(t) = a*t^(n-1)*cos(2*pi*fc*t)*exp(-2*pi*beta*t)
%
% .. math:: g(t) = at^{n-1}cos(2\pi\cdot fc\cdot t)e^{-2\pi \beta \cdot t}
%
% The gammatone filters as implemented by this function generate
% complex valued output, because the filters are modulated by the
% exponential function. Using `real` on the output will give the
% coefficients of the corresponding cosine modulated filters.
%
% To create the filter coefficients of a 1-erb spaced filter bank using
% gammatone filters use the following construction::
%
% g = gammatonefir(erbspacebw(flow,fhigh),fs);
%
% See also: erbspace, audspace, audfiltbw
%
% Demos: demo_auditoryfilterbank
%
% References: aertsen1980strI glasberg1990daf
% AUTHOR : Peter L. S��ndergaard
% ------ Checking of input parameters ---------
if nargin<2
error('Too few input arguments.');
end;
if ~isnumeric(fs) || ~isscalar(fs) || fs<=0
error('%s: fs must be a positive scalar.',upper(mfilename));
end;
if ~isnumeric(fc) || ~isvector(fc) || any(fc<0) || any(fc>fs/2)
error(['%s: fc must be a vector of positive values that are less than half ' ...
'the sampling rate.'],upper(mfilename));
end;
definput.import={'normalize'};
definput.importdefaults={'null'};
definput.flags.real={'complex','real'};
definput.keyvals.n=[];
definput.flags.phase={'causalphase','peakphase'};
definput.keyvals.betamul=1.0183;
[flags,keyvals,n,betamul] = ltfatarghelper({'n','betamul'},definput,varargin);
nchannels = length(fc);
% ourbeta is used in order not to mask the beta function.
ourbeta = betamul*audfiltbw(fc);
if isempty(n)
% Calculate a good value for n
% FIXME actually do this
n=5000;
end;
b={};
for ii = 1:nchannels
delay = 3/(2*pi*ourbeta(ii));
scalconst = 2*(2*pi*ourbeta(ii))^4/factorial(4-1)/fs;
nfirst = ceil(fs*delay);
if nfirst>n/2
error(['%s: The desired filter length is too short to accomodate the ' ...
'beginning of the filter. Please choose a filter length of ' ...
'at least %i samples.'],upper(mfilename),nfirst*2);
end;
nlast = n/2;
t=[(0:nlast-1)/fs+delay,...
(0:nfirst-1)/fs-nfirst/fs+delay].';
% g(t) = a*t^(n-1)*cos(2*pi*fc*t)*exp(-2*pi*beta*t)
if flags.do_real
bwork = scalconst*t.^(4-1).*cos(2*pi*fc(ii)*t).*exp(-2*pi* ...
ourbeta(ii)*t);
else
bwork = scalconst*t.^(4-1).*exp(2*pi*i*fc(ii)*t).*exp(-2*pi* ...
ourbeta(ii)*t);
end;
if flags.do_peakphase
bwork=bwork*exp(-2*pi*i*fc(ii)*delay);
end;
% Insert zeros before the start of the signal.
b{ii}=[bwork(1:nlast);zeros(n-nlast-nfirst,1);bwork(nlast+1:nlast+nfirst)];
b{ii}=normalize(b{ii},flags.norm);
end;