Mercurial > hg > camir-aes2014

annotate toolboxes/MIRtoolbox1.3.2/MIRToolbox/mirfilterbank.m @ 0:e9a9cd732c1e tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
first hg version after svn
author wolffd
date Tue, 10 Feb 2015 15:05:51 +0000
parents
children
rev   line source
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}]}});