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
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
|
function b=gammatonefir(fc,fs,varargin);
%-*- texinfo -*-
%@deftypefn {Function} gammatonefir
%@verbatim
%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 : max. filter length.
% beta : bandwidth of the filter.
%
% Output parameters:
% b : FIR filters as an cell-array of structs.
%
% GAMMATONEFIR(fc,fs,n,betamul) computes the filter coefficients of a
% digital FIR gammatone filter with length at most n, 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. The actual length of the inpulse response
% depends on fc (the filter is longer for low center frequencies),
% fs and betamul but it is never bigger than n.
%
% 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^(4-1)*cos(2*pi*fc*t)*exp(-2*pi*beta*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);
%
%
%
% References:
% A. Aertsen and P. Johannesma. Spectro-temporal receptive fields of
% auditory neurons in the grassfrog. I. Characterization of tonal and
% natural stimuli. Biol. Cybern, 38:223--234, 1980.
%
% B. R. Glasberg and B. Moore. Derivation of auditory filter shapes from
% notched-noise data. Hearing Research, 47(1-2):103, 1990.
%
%@end verbatim
%@strong{Url}: @url{http://ltfat.github.io/doc/auditory/gammatonefir.html}
%@seealso{erbspace, audspace, audfiltbw, demo_auditoryfilterbank}
%@end deftypefn
% Copyright (C) 2005-2016 Peter L. Soendergaard <peter@sonderport.dk>.
% This file is part of LTFAT version 2.3.1
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
% AUTHOR : Peter L. Soendergaard
% ------ 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=cell(nchannels,1);
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 = floor(n/2);
t=[(0:nfirst-1)/fs-nfirst/fs+delay,(0:nlast-1)/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.
%bwork = fftshift([bwork(1:nlast);zeros(n-nlast-nfirst,1);bwork(nlast+1:nlast+nfirst)]);
bwork = normalize(bwork,flags.norm);
b{ii} = struct('h',bwork,'offset',-nfirst,'realonly',0);
end;
|
