wolffd@0
|
1 function varargout = mirfilterbank(orig,varargin)
|
wolffd@0
|
2 % b = mirfilterbank(x) performs a filterbank decomposition of an audio
|
wolffd@0
|
3 % waveform.
|
wolffd@0
|
4 % Optional arguments:
|
wolffd@0
|
5 % mirfilterbank(...,t) selects a type of filterbank.
|
wolffd@0
|
6 % Possible values:
|
wolffd@0
|
7 % t = 'Gammatone' (default for audio files).
|
wolffd@0
|
8 % Requires the Auditory Toolbox.
|
wolffd@0
|
9 % mirfilterbank(...,'Lowest',f): lowest frequency in Hz
|
wolffd@0
|
10 % (default: 50)
|
wolffd@0
|
11 % t = '2Channels' proposed in (Tolonen & Karjalainen, 2000)
|
wolffd@0
|
12 % mirfilterbank(...,'NbChannels',N), or simply filterbank(x,N):
|
wolffd@0
|
13 % specifies the number of channels in the bank.
|
wolffd@0
|
14 % (default: N = 10)
|
wolffd@0
|
15 % mirfilterbank(...,'Channel',c) only output the channels whose
|
wolffd@0
|
16 % ranks are indicated in the array c.
|
wolffd@0
|
17 % (default: c = (1:N))
|
wolffd@0
|
18 % mirfilterbank(...,'Manual',f) specifies a set of low-pass, band-
|
wolffd@0
|
19 % pass and high-pass eliptic filters (Scheirer, 1998).
|
wolffd@0
|
20 % The series of cut-off frequencies as to be specified
|
wolffd@0
|
21 % as next parameter f.
|
wolffd@0
|
22 % If this series of frequencies f begins with -Inf,
|
wolffd@0
|
23 % the first filter is low-pass.
|
wolffd@0
|
24 % If this series of frequencies f ends with Inf,
|
wolffd@0
|
25 % the last filter is high-pass.
|
wolffd@0
|
26 % mirfilterbank(...,'Order',o) specifies the order of the filter.
|
wolffd@0
|
27 % (Default: o = 4) (Scheirer, 1998)
|
wolffd@0
|
28 % mirfilterbank(...,'Hop',h) specifies the degree of spectral
|
wolffd@0
|
29 % overlapping between successive channels.
|
wolffd@0
|
30 % If h = 1 (default value), the filters are non-overlapping.
|
wolffd@0
|
31 % If h = 2, the filters are half-overlapping.
|
wolffd@0
|
32 % If h = 3, the spectral hop factor between successive
|
wolffd@0
|
33 % filters is a third of the whole frequency region, etc.
|
wolffd@0
|
34 % mirfilterbank(...,p) specifies predefined filterbanks, all
|
wolffd@0
|
35 % implemented using elliptic filters, by default of order 4.
|
wolffd@0
|
36 % Possible values:
|
wolffd@0
|
37 % p = 'Mel' (mel-scale)
|
wolffd@0
|
38 % p = 'Bark' (bark-scale)
|
wolffd@0
|
39 % p = 'Scheirer' proposed in (Scheirer, 1998) corresponds to
|
wolffd@0
|
40 % 'Manual',[-Inf 200 400 800 1600 3200 Inf]
|
wolffd@0
|
41 % p = 'Klapuri' proposed in (Klapuri, 1999) corresponds to
|
wolffd@0
|
42 % 'Manual',44*[2.^ ([ 0:2, ( 9+(0:17) )/3 ]) ]
|
wolffd@0
|
43
|
wolffd@0
|
44 nCh.key = 'NbChannels';
|
wolffd@0
|
45 nCh.type = 'Integer';
|
wolffd@0
|
46 nCh.default = 10;
|
wolffd@0
|
47 option.nCh = nCh;
|
wolffd@0
|
48
|
wolffd@0
|
49 Ch.key = {'Channel','Channels'};
|
wolffd@0
|
50 Ch.type = 'Integer';
|
wolffd@0
|
51 Ch.default = 0;
|
wolffd@0
|
52 option.Ch = Ch;
|
wolffd@0
|
53
|
wolffd@0
|
54 lowF.key = 'Lowest';
|
wolffd@0
|
55 lowF.type = 'Integer';
|
wolffd@0
|
56 lowF.default = 50;
|
wolffd@0
|
57 option.lowF = lowF;
|
wolffd@0
|
58
|
wolffd@0
|
59 freq.key = 'Manual';
|
wolffd@0
|
60 freq.type = 'Integer';
|
wolffd@0
|
61 freq.default = NaN;
|
wolffd@0
|
62 option.freq = freq;
|
wolffd@0
|
63
|
wolffd@0
|
64 overlap.key = 'Hop';
|
wolffd@0
|
65 overlap.type = 'Boolean';
|
wolffd@0
|
66 overlap.default = 1;
|
wolffd@0
|
67 option.overlap = overlap;
|
wolffd@0
|
68
|
wolffd@0
|
69 filtertype.type = 'String';
|
wolffd@0
|
70 filtertype.choice = {'Gammatone','2Channels',0};
|
wolffd@0
|
71 filtertype.default = 'Gammatone';
|
wolffd@0
|
72 option.filtertype = filtertype;
|
wolffd@0
|
73
|
wolffd@0
|
74 filterorder.key = 'Order';
|
wolffd@0
|
75 filterorder.type = 'Integer';
|
wolffd@0
|
76 filterorder.default = 4;
|
wolffd@0
|
77 option.filterorder = filterorder;
|
wolffd@0
|
78
|
wolffd@0
|
79 presel.type = 'String';
|
wolffd@0
|
80 presel.choice = {'Scheirer','Klapuri','Mel','Bark'};
|
wolffd@0
|
81 presel.default = '';
|
wolffd@0
|
82 option.presel = presel;
|
wolffd@0
|
83
|
wolffd@0
|
84 specif.option = option;
|
wolffd@0
|
85
|
wolffd@0
|
86 specif.eachchunk = @eachchunk;
|
wolffd@0
|
87 specif.combinechunk = @combinechunk;
|
wolffd@0
|
88
|
wolffd@0
|
89 varargout = mirfunction(@mirfilterbank,orig,varargin,nargout,specif,@init,@main);
|
wolffd@0
|
90
|
wolffd@0
|
91
|
wolffd@0
|
92 function [x type] = init(x,option)
|
wolffd@0
|
93 type = 'miraudio';
|
wolffd@0
|
94
|
wolffd@0
|
95
|
wolffd@0
|
96 function b = main(x,option,postoption)
|
wolffd@0
|
97 if iscell(x)
|
wolffd@0
|
98 x = x{1};
|
wolffd@0
|
99 end
|
wolffd@0
|
100 f = get(x,'Sampling');
|
wolffd@0
|
101
|
wolffd@0
|
102 if strcmpi(option.presel,'Scheirer')
|
wolffd@0
|
103 option.freq = [-Inf 200 400 800 1600 3200 Inf];
|
wolffd@0
|
104 elseif strcmpi(option.presel,'Klapuri')
|
wolffd@0
|
105 option.freq = 44*[2.^([0:2,(9+(0:17))/3])];
|
wolffd@0
|
106 elseif strcmpi(option.presel,'Mel')
|
wolffd@0
|
107 lowestFrequency = 133.3333;
|
wolffd@0
|
108 linearFilters = 13;
|
wolffd@0
|
109 linearSpacing = 66.66666666;
|
wolffd@0
|
110 logFilters = 27;
|
wolffd@0
|
111 logSpacing = 1.0711703;
|
wolffd@0
|
112 totalFilters = linearFilters + logFilters;
|
wolffd@0
|
113 cepstralCoefficients = 13;
|
wolffd@0
|
114 option.freq = lowestFrequency + (0:linearFilters-1)*linearSpacing;
|
wolffd@0
|
115 option.freq(linearFilters+1:totalFilters+2) = ...
|
wolffd@0
|
116 option.freq(linearFilters) * logSpacing.^(1:logFilters+2);
|
wolffd@0
|
117
|
wolffd@0
|
118 option.overlap = 2;
|
wolffd@0
|
119 elseif strcmpi(option.presel,'Bark')
|
wolffd@0
|
120 option.freq = [10 20 30 40 51 63 77 92 108 127 148 172 200 232 ...
|
wolffd@0
|
121 270 315 370 440 530 640 770 950 1200 1550]*10; %% Hz
|
wolffd@0
|
122 end
|
wolffd@0
|
123 if not(isnan(option.freq))
|
wolffd@0
|
124 option.filtertype = 'Manual';
|
wolffd@0
|
125 end
|
wolffd@0
|
126
|
wolffd@0
|
127 d = get(x,'Data');
|
wolffd@0
|
128 if size(d{1}{1},3) > 1
|
wolffd@0
|
129 warning('WARNING IN MIRFILTERBANK: The input data is already decomposed into channels. No more channel decomposition.');
|
wolffd@0
|
130 if option.Ch
|
wolffd@0
|
131 cx = get(x,'Channels');
|
wolffd@0
|
132 db = cell(1,length(d));
|
wolffd@0
|
133 for k = 1:length(d)
|
wolffd@0
|
134 db{k} = cell(1,length(d{k}));
|
wolffd@0
|
135 for l = 1:length(d{k})
|
wolffd@0
|
136 for i = 1:length(option.Ch)
|
wolffd@0
|
137 db{k}{l}(:,:,i) = d{k}{l}(:,:,find(cx{k} == option.Ch(i)));
|
wolffd@0
|
138 end
|
wolffd@0
|
139 end
|
wolffd@0
|
140 end
|
wolffd@0
|
141 b = set(x,'Data',db,'Channels',option.Ch);
|
wolffd@0
|
142 else
|
wolffd@0
|
143 b = x;
|
wolffd@0
|
144 end
|
wolffd@0
|
145 else
|
wolffd@0
|
146 i = 1;
|
wolffd@0
|
147 while i <= length(d)
|
wolffd@0
|
148 if size(d{i}{1},2) > 1
|
wolffd@0
|
149 warning('WARNING IN MIRFILTERBANK: The frame decomposition should be performed after the filterbank decomposition.');
|
wolffd@0
|
150 disp(['Suggestion: Use the ' char(34) 'Frame' char(34) ' option available to some of the MIRtoolbox functions.'])
|
wolffd@0
|
151 break
|
wolffd@0
|
152 end
|
wolffd@0
|
153 i = i+1;
|
wolffd@0
|
154 end
|
wolffd@0
|
155 nCh = option.nCh;
|
wolffd@0
|
156 Ch = option.Ch;
|
wolffd@0
|
157 if strcmpi(option.filtertype, 'Gammatone');
|
wolffd@0
|
158 [tmp x] = gettmp(x); %get(x,'Tmp');
|
wolffd@0
|
159 if not(Ch)
|
wolffd@0
|
160 Ch = 1:nCh;
|
wolffd@0
|
161 end
|
wolffd@0
|
162 erb = cell(1,length(d));
|
wolffd@0
|
163 for i = 1:length(d)
|
wolffd@0
|
164 erb{i} = cell(1,length(d{i}));
|
wolffd@0
|
165 nch{i} = Ch;
|
wolffd@0
|
166 for j = 1:length(d{i})
|
wolffd@0
|
167 try
|
wolffd@0
|
168 coef = MakeERBFilters(f{i},nCh,option.lowF);
|
wolffd@0
|
169 catch
|
wolffd@0
|
170 error(['ERROR IN MIRFILTERBANK: Auditory Toolbox needs to be installed.']);
|
wolffd@0
|
171 end
|
wolffd@0
|
172 [erb{i}{j} tmp] = ERBfilterbank(d{i}{j},coef,Ch,tmp);
|
wolffd@0
|
173 end
|
wolffd@0
|
174 end
|
wolffd@0
|
175 b = set(x,'Data',erb,'Channels',nch);%,'Tmp',tmp);
|
wolffd@0
|
176 clear erb
|
wolffd@0
|
177 b = settmp(b,tmp);
|
wolffd@0
|
178 elseif strcmpi(option.filtertype, '2Channels');
|
wolffd@0
|
179 if not(Ch)
|
wolffd@0
|
180 Ch = 1:2;
|
wolffd@0
|
181 end
|
wolffd@0
|
182 output = cell(1,length(d));
|
wolffd@0
|
183 try
|
wolffd@0
|
184 for i = 1:length(d)
|
wolffd@0
|
185 output{i} = cell(1,length(d{i}));
|
wolffd@0
|
186 [bl,al] = butter(4,[70 1000]/f{i}*2);
|
wolffd@0
|
187 k = 1;
|
wolffd@0
|
188 if ismember(1,Ch)
|
wolffd@0
|
189 Hdl = dfilt.df2t(bl,al);
|
wolffd@0
|
190 %fvtool(Hdl)
|
wolffd@0
|
191 Hdl.PersistentMemory = true;
|
wolffd@0
|
192 for j = 1:length(d{i})
|
wolffd@0
|
193 low = filter(Hdl,d{i}{j});
|
wolffd@0
|
194 output{i}{j}(:,:,k) = low;
|
wolffd@0
|
195 end
|
wolffd@0
|
196 k = k+1;
|
wolffd@0
|
197 end
|
wolffd@0
|
198 if ismember(2,Ch)
|
wolffd@0
|
199 if f{i} < 20000
|
wolffd@0
|
200 [bh,ah] = butter(2,1000/f{i}*2,'high');
|
wolffd@0
|
201 else
|
wolffd@0
|
202 [bh,ah] = butter(2,[1000 10000]/f{i}*2);
|
wolffd@0
|
203 end
|
wolffd@0
|
204 Hdlh = dfilt.df2t(bl,al);
|
wolffd@0
|
205 %fvtool(Hdlh)
|
wolffd@0
|
206 Hdh = dfilt.df2t(bh,ah);
|
wolffd@0
|
207 %fvtool(Hdh)
|
wolffd@0
|
208 Hdlh.PersistentMemory = true;
|
wolffd@0
|
209 Hdh.PersistentMemory = true;
|
wolffd@0
|
210 for j = 1:length(d{i})
|
wolffd@0
|
211 high = filter(Hdh,d{i}{j});
|
wolffd@0
|
212 high = max(high,0);
|
wolffd@0
|
213 high = filter(Hdlh,high);
|
wolffd@0
|
214 output{i}{j}(:,:,k) = high;
|
wolffd@0
|
215 end
|
wolffd@0
|
216 end
|
wolffd@0
|
217 nch{i} = Ch;
|
wolffd@0
|
218 end
|
wolffd@0
|
219 b = set(x,'Data',output,'Channels',nch);
|
wolffd@0
|
220 catch
|
wolffd@0
|
221 warning(['WARNING IN MIRFILTERBANK: Signal Processing Toolbox (version 6.2.1, or higher) not installed: no filterbank decomposition.']);
|
wolffd@0
|
222 b = x;
|
wolffd@0
|
223 end
|
wolffd@0
|
224 elseif strcmpi(option.filtertype,'Manual');
|
wolffd@0
|
225 output = cell(1,length(d));
|
wolffd@0
|
226 for i = 1:length(d)
|
wolffd@0
|
227 freqi = option.freq;
|
wolffd@0
|
228 j = 1;
|
wolffd@0
|
229 while j <= length(freqi)
|
wolffd@0
|
230 if not(isinf(freqi(j))) && freqi(j)>f{i}/2
|
wolffd@0
|
231 if j == length(freqi)
|
wolffd@0
|
232 freqi(j) = Inf;
|
wolffd@0
|
233 else
|
wolffd@0
|
234 freqi(j) = [];
|
wolffd@0
|
235 j = j-1;
|
wolffd@0
|
236 end
|
wolffd@0
|
237 end
|
wolffd@0
|
238 j = j+1;
|
wolffd@0
|
239 end
|
wolffd@0
|
240 output{i} = cell(1,length(d{i}));
|
wolffd@0
|
241 step = option.overlap;
|
wolffd@0
|
242 for j = 1:length(freqi)-step
|
wolffd@0
|
243 if isinf(freqi(j))
|
wolffd@0
|
244 [z{j},p{j},k{j}] = ellip(option.filterorder,3,40,...
|
wolffd@0
|
245 freqi(j+step)/f{i}*2);
|
wolffd@0
|
246 elseif isinf(freqi(j+step))
|
wolffd@0
|
247 [z{j},p{j},k{j}] = ellip(option.filterorder,3,40,...
|
wolffd@0
|
248 freqi(j)/f{i}*2,'high');
|
wolffd@0
|
249 else
|
wolffd@0
|
250 [z{j},p{j},k{j}] = ellip(option.filterorder,3,40,...
|
wolffd@0
|
251 freqi([j j+step])/f{i}*2);
|
wolffd@0
|
252 end
|
wolffd@0
|
253 end
|
wolffd@0
|
254 for j = 1:length(z)
|
wolffd@0
|
255 [sos,g] = zp2sos(z{j},p{j},k{j});
|
wolffd@0
|
256 Hd(j) = dfilt.df2tsos(sos,g);
|
wolffd@0
|
257 Hd(j).PersistentMemory = true;
|
wolffd@0
|
258 for h = 1:length(d{i})
|
wolffd@0
|
259 output{i}{h}(:,:,j) = filter(Hd(j),d{i}{h});
|
wolffd@0
|
260 end
|
wolffd@0
|
261 end
|
wolffd@0
|
262 %fvtool(Hd)
|
wolffd@0
|
263 nch{i} = 1:length(freqi)-step;
|
wolffd@0
|
264 end
|
wolffd@0
|
265 b = set(x,'Data',output,'Channels',nch);
|
wolffd@0
|
266 end
|
wolffd@0
|
267 end
|
wolffd@0
|
268
|
wolffd@0
|
269
|
wolffd@0
|
270 function [output tmp] = ERBfilterbank(x, fcoefs, chans, tmp)
|
wolffd@0
|
271 % The following is based on the source code from Auditory Toolbox
|
wolffd@0
|
272 % (A part that I could not call directly from MIRtoolbox, and that I
|
wolffd@0
|
273 % generalize to chunk decomposition)
|
wolffd@0
|
274
|
wolffd@0
|
275 % (Malcolm Slaney, August 1993, (c) 1998 Interval Research Corporation)
|
wolffd@0
|
276
|
wolffd@0
|
277 % Process an input waveform with a gammatone filter bank.
|
wolffd@0
|
278 % The fcoefs parameter, which completely specifies the Gammatone filterbank,
|
wolffd@0
|
279 % should be designed with the MakeERBFilters function.
|
wolffd@0
|
280
|
wolffd@0
|
281 A0 = fcoefs(:,1);
|
wolffd@0
|
282 A11 = fcoefs(:,2);
|
wolffd@0
|
283 A12 = fcoefs(:,3);
|
wolffd@0
|
284 A13 = fcoefs(:,4);
|
wolffd@0
|
285 A14 = fcoefs(:,5);
|
wolffd@0
|
286 A2 = fcoefs(:,6);
|
wolffd@0
|
287 B0 = fcoefs(:,7);
|
wolffd@0
|
288 B1 = fcoefs(:,8);
|
wolffd@0
|
289 B2 = fcoefs(:,9);
|
wolffd@0
|
290 gain= fcoefs(:,10);
|
wolffd@0
|
291 lc = length(chans);
|
wolffd@0
|
292 output = zeros(size(x,1),size(x,2),lc);
|
wolffd@0
|
293 if isempty(tmp)
|
wolffd@0
|
294 emptytmp = 1;
|
wolffd@0
|
295 else
|
wolffd@0
|
296 emptytmp = 0;
|
wolffd@0
|
297 end
|
wolffd@0
|
298 for i = 1:lc
|
wolffd@0
|
299 chan = length(gain)-chans(i)+1; % Revert the channels order
|
wolffd@0
|
300 aa1 = [A0(chan)/gain(chan) A11(chan)/gain(chan) A2(chan)/gain(chan)];
|
wolffd@0
|
301 bb1 = [B0(chan) B1(chan) B2(chan)];
|
wolffd@0
|
302 aa2 = [A0(chan) A12(chan) A2(chan)];
|
wolffd@0
|
303 bb2 = [B0(chan) B1(chan) B2(chan)];
|
wolffd@0
|
304 aa3 = [A0(chan) A12(chan) A2(chan)];
|
wolffd@0
|
305 bb3 = [B0(chan) B1(chan) B2(chan)];
|
wolffd@0
|
306 aa4 = [A0(chan) A13(chan) A2(chan)];
|
wolffd@0
|
307 bb4 = [B0(chan) B1(chan) B2(chan)];
|
wolffd@0
|
308 if emptytmp
|
wolffd@0
|
309 [y1 tmp(:,:,i,1)] = filter(aa1,bb1,x);
|
wolffd@0
|
310 [y2 tmp(:,:,i,2)] = filter(aa2,bb2,y1);
|
wolffd@0
|
311 [y3 tmp(:,:,i,3)] = filter(aa3,bb3,y2);
|
wolffd@0
|
312 [y4 tmp(:,:,i,4)] = filter(aa4,bb4,y3);
|
wolffd@0
|
313 else
|
wolffd@0
|
314 [y1 tmp(:,:,i,1)] = filter(aa1,bb1,x,tmp(:,:,i,1));
|
wolffd@0
|
315 [y2 tmp(:,:,i,2)] = filter(aa2,bb2,y1,tmp(:,:,i,2));
|
wolffd@0
|
316 [y3 tmp(:,:,i,3)] = filter(aa3,bb3,y2,tmp(:,:,i,3));
|
wolffd@0
|
317 [y4 tmp(:,:,i,4)] = filter(aa4,bb4,y3,tmp(:,:,i,4));
|
wolffd@0
|
318 end
|
wolffd@0
|
319 output(:,:,i) = y4;
|
wolffd@0
|
320 end
|
wolffd@0
|
321
|
wolffd@0
|
322
|
wolffd@0
|
323 function [y orig] = eachchunk(orig,option,missing)
|
wolffd@0
|
324 y = mirfilterbank(orig,option);
|
wolffd@0
|
325
|
wolffd@0
|
326
|
wolffd@0
|
327 function y = combinechunk(old,new)
|
wolffd@0
|
328 do = get(old,'Data');
|
wolffd@0
|
329 to = get(old,'Time');
|
wolffd@0
|
330 dn = get(new,'Data');
|
wolffd@0
|
331 tn = get(new,'Time');
|
wolffd@0
|
332 y = set(old,'Data',{{[do{1}{1};dn{1}{1}]}},...
|
wolffd@0
|
333 'Time',{{[to{1}{1};tn{1}{1}]}}); |