wolffd@0: function varargout = mirfilterbank(orig,varargin)
wolffd@0: %   b = mirfilterbank(x) performs a filterbank decomposition of an audio
wolffd@0: %       waveform.
wolffd@0: %   Optional arguments:
wolffd@0: %       mirfilterbank(...,t) selects a type of filterbank.
wolffd@0: %           Possible values:
wolffd@0: %               t = 'Gammatone' (default for audio files).
wolffd@0: %                   Requires the Auditory Toolbox.
wolffd@0: %                   mirfilterbank(...,'Lowest',f): lowest frequency in Hz
wolffd@0: %                       (default: 50)
wolffd@0: %               t = '2Channels' proposed in (Tolonen & Karjalainen, 2000)
wolffd@0: %           mirfilterbank(...,'NbChannels',N), or simply filterbank(x,N):
wolffd@0: %               specifies the number of channels in the bank.
wolffd@0: %               (default: N = 10)
wolffd@0: %           mirfilterbank(...,'Channel',c) only output the channels whose
wolffd@0: %               ranks are indicated in the array c.
wolffd@0: %               (default: c = (1:N))
wolffd@0: %       mirfilterbank(...,'Manual',f) specifies a set of low-pass, band-
wolffd@0: %                   pass and high-pass eliptic filters (Scheirer, 1998). 
wolffd@0: %                   The series of cut-off frequencies as to be specified
wolffd@0: %                   as next parameter f.
wolffd@0: %                       If this series of frequencies f begins with -Inf,
wolffd@0: %                           the first filter is low-pass.
wolffd@0: %                       If this series of frequencies f ends with Inf,
wolffd@0: %                           the last filter is high-pass.
wolffd@0: %           mirfilterbank(...,'Order',o) specifies the order of the filter.
wolffd@0: %               (Default: o = 4) (Scheirer, 1998)
wolffd@0: %           mirfilterbank(...,'Hop',h) specifies the degree of spectral
wolffd@0: %               overlapping between successive channels.
wolffd@0: %               If h = 1 (default value), the filters are non-overlapping.
wolffd@0: %               If h = 2, the filters are half-overlapping.
wolffd@0: %               If h = 3, the spectral hop factor between successive
wolffd@0: %                   filters is a third of the whole frequency region, etc.
wolffd@0: %       mirfilterbank(...,p) specifies predefined filterbanks, all
wolffd@0: %           implemented using elliptic filters, by default of order 4.
wolffd@0: %           Possible values:
wolffd@0: %               p = 'Mel' (mel-scale)
wolffd@0: %               p = 'Bark' (bark-scale)
wolffd@0: %               p = 'Scheirer' proposed in (Scheirer, 1998) corresponds to
wolffd@0: %                   'Manual',[-Inf 200 400 800 1600 3200 Inf]
wolffd@0: %               p = 'Klapuri' proposed in (Klapuri, 1999) corresponds to
wolffd@0: %                   'Manual',44*[2.^ ([ 0:2, ( 9+(0:17) )/3 ]) ]
wolffd@0: 
wolffd@0:         nCh.key = 'NbChannels';
wolffd@0:         nCh.type = 'Integer';
wolffd@0:         nCh.default = 10;
wolffd@0:     option.nCh = nCh;
wolffd@0:     
wolffd@0:         Ch.key = {'Channel','Channels'};
wolffd@0:         Ch.type = 'Integer';
wolffd@0:         Ch.default = 0;
wolffd@0:     option.Ch = Ch;
wolffd@0:     
wolffd@0:         lowF.key = 'Lowest';
wolffd@0:         lowF.type = 'Integer';
wolffd@0:         lowF.default = 50;
wolffd@0:     option.lowF = lowF;
wolffd@0: 
wolffd@0:         freq.key = 'Manual';
wolffd@0:         freq.type = 'Integer';
wolffd@0:         freq.default = NaN;
wolffd@0:     option.freq = freq;
wolffd@0:     
wolffd@0:         overlap.key = 'Hop';
wolffd@0:         overlap.type = 'Boolean';
wolffd@0:         overlap.default = 1;
wolffd@0:     option.overlap = overlap;
wolffd@0: 
wolffd@0:         filtertype.type = 'String';
wolffd@0:         filtertype.choice = {'Gammatone','2Channels',0};
wolffd@0:         filtertype.default = 'Gammatone';
wolffd@0:     option.filtertype = filtertype;
wolffd@0: 
wolffd@0:         filterorder.key = 'Order';
wolffd@0:         filterorder.type = 'Integer';
wolffd@0:         filterorder.default = 4;
wolffd@0:     option.filterorder = filterorder;
wolffd@0: 
wolffd@0:         presel.type = 'String';
wolffd@0:         presel.choice = {'Scheirer','Klapuri','Mel','Bark'};
wolffd@0:         presel.default = '';
wolffd@0:     option.presel = presel;
wolffd@0: 
wolffd@0: specif.option = option;
wolffd@0: 
wolffd@0: specif.eachchunk = @eachchunk;
wolffd@0: specif.combinechunk = @combinechunk;
wolffd@0: 
wolffd@0: varargout = mirfunction(@mirfilterbank,orig,varargin,nargout,specif,@init,@main);
wolffd@0: 
wolffd@0: 
wolffd@0: function [x type] = init(x,option)
wolffd@0: type = 'miraudio';
wolffd@0: 
wolffd@0: 
wolffd@0: function b = main(x,option,postoption)
wolffd@0: if iscell(x)
wolffd@0:     x = x{1};
wolffd@0: end
wolffd@0: f = get(x,'Sampling');
wolffd@0:     
wolffd@0: if strcmpi(option.presel,'Scheirer')
wolffd@0:     option.freq = [-Inf 200 400 800 1600 3200 Inf];
wolffd@0: elseif strcmpi(option.presel,'Klapuri')
wolffd@0:     option.freq = 44*[2.^([0:2,(9+(0:17))/3])];
wolffd@0: elseif strcmpi(option.presel,'Mel')
wolffd@0:     lowestFrequency = 133.3333;
wolffd@0:     linearFilters = 13;
wolffd@0:     linearSpacing = 66.66666666;
wolffd@0:     logFilters = 27;
wolffd@0:     logSpacing = 1.0711703;
wolffd@0:     totalFilters = linearFilters + logFilters;
wolffd@0:     cepstralCoefficients = 13;
wolffd@0:     option.freq = lowestFrequency + (0:linearFilters-1)*linearSpacing;
wolffd@0:     option.freq(linearFilters+1:totalFilters+2) = ...
wolffd@0:         option.freq(linearFilters) * logSpacing.^(1:logFilters+2);
wolffd@0: 
wolffd@0:     option.overlap = 2;
wolffd@0: elseif strcmpi(option.presel,'Bark')
wolffd@0:     option.freq = [10 20 30 40 51 63 77 92 108 127 148 172 200 232 ...
wolffd@0:                     270 315 370 440 530 640 770 950 1200 1550]*10; %% Hz
wolffd@0: end
wolffd@0: if not(isnan(option.freq))
wolffd@0:     option.filtertype = 'Manual';
wolffd@0: end
wolffd@0:         
wolffd@0: d = get(x,'Data');
wolffd@0: if size(d{1}{1},3) > 1
wolffd@0:     warning('WARNING IN MIRFILTERBANK: The input data is already decomposed into channels. No more channel decomposition.');
wolffd@0:     if option.Ch
wolffd@0:         cx = get(x,'Channels');
wolffd@0:         db = cell(1,length(d));
wolffd@0:         for k = 1:length(d)
wolffd@0:             db{k} = cell(1,length(d{k}));
wolffd@0:             for l = 1:length(d{k})
wolffd@0:                 for i = 1:length(option.Ch)
wolffd@0:                     db{k}{l}(:,:,i) = d{k}{l}(:,:,find(cx{k} == option.Ch(i)));
wolffd@0:                 end
wolffd@0:             end
wolffd@0:         end
wolffd@0:         b = set(x,'Data',db,'Channels',option.Ch);
wolffd@0:     else
wolffd@0:         b = x;
wolffd@0:     end
wolffd@0: else
wolffd@0:     i = 1;
wolffd@0:     while i <= length(d)
wolffd@0:         if size(d{i}{1},2) > 1
wolffd@0:             warning('WARNING IN MIRFILTERBANK: The frame decomposition should be performed after the filterbank decomposition.');
wolffd@0:             disp(['Suggestion: Use the ' char(34) 'Frame' char(34) ' option available to some of the MIRtoolbox functions.'])
wolffd@0:             break
wolffd@0:         end
wolffd@0:         i = i+1;
wolffd@0:     end
wolffd@0:     nCh = option.nCh;
wolffd@0:     Ch = option.Ch;
wolffd@0:     if strcmpi(option.filtertype, 'Gammatone');
wolffd@0:         [tmp x] = gettmp(x); %get(x,'Tmp');
wolffd@0:         if not(Ch)
wolffd@0:             Ch = 1:nCh;
wolffd@0:         end
wolffd@0:         erb = cell(1,length(d));
wolffd@0:         for i = 1:length(d)
wolffd@0:             erb{i} = cell(1,length(d{i}));
wolffd@0:             nch{i} = Ch;
wolffd@0:             for j = 1:length(d{i})
wolffd@0:                 try
wolffd@0:                     coef = MakeERBFilters(f{i},nCh,option.lowF);
wolffd@0:                 catch
wolffd@0:                     error(['ERROR IN MIRFILTERBANK: Auditory Toolbox needs to be installed.']);
wolffd@0:                 end                    
wolffd@0:                 [erb{i}{j} tmp] = ERBfilterbank(d{i}{j},coef,Ch,tmp);
wolffd@0:             end
wolffd@0:         end
wolffd@0:         b = set(x,'Data',erb,'Channels',nch);%,'Tmp',tmp);
wolffd@0:         clear erb
wolffd@0:         b = settmp(b,tmp);
wolffd@0:     elseif strcmpi(option.filtertype, '2Channels');
wolffd@0:         if not(Ch)
wolffd@0:             Ch = 1:2;
wolffd@0:         end
wolffd@0:         output = cell(1,length(d));
wolffd@0:         try
wolffd@0:             for i = 1:length(d)
wolffd@0:                 output{i} = cell(1,length(d{i}));
wolffd@0:                 [bl,al] = butter(4,[70 1000]/f{i}*2);
wolffd@0:                 k = 1;
wolffd@0:                 if ismember(1,Ch)
wolffd@0:                     Hdl = dfilt.df2t(bl,al);
wolffd@0:                     %fvtool(Hdl)
wolffd@0:                     Hdl.PersistentMemory = true;
wolffd@0:                     for j = 1:length(d{i})
wolffd@0:                         low = filter(Hdl,d{i}{j});
wolffd@0:                         output{i}{j}(:,:,k) = low;
wolffd@0:                     end
wolffd@0:                     k = k+1;
wolffd@0:                 end
wolffd@0:                 if ismember(2,Ch)
wolffd@0:                     if f{i} < 20000
wolffd@0:                         [bh,ah] = butter(2,1000/f{i}*2,'high');
wolffd@0:                     else
wolffd@0:                         [bh,ah] = butter(2,[1000 10000]/f{i}*2);
wolffd@0:                     end
wolffd@0:                     Hdlh = dfilt.df2t(bl,al);
wolffd@0:                     %fvtool(Hdlh)
wolffd@0:                     Hdh = dfilt.df2t(bh,ah);
wolffd@0:                     %fvtool(Hdh)
wolffd@0:                     Hdlh.PersistentMemory = true;
wolffd@0:                     Hdh.PersistentMemory = true;
wolffd@0:                     for j = 1:length(d{i})
wolffd@0:                         high = filter(Hdh,d{i}{j});
wolffd@0:                         high = max(high,0);
wolffd@0:                         high = filter(Hdlh,high);
wolffd@0:                         output{i}{j}(:,:,k) = high;
wolffd@0:                     end
wolffd@0:                 end
wolffd@0:                 nch{i} = Ch;
wolffd@0:             end
wolffd@0:             b = set(x,'Data',output,'Channels',nch);
wolffd@0:         catch
wolffd@0:             warning(['WARNING IN MIRFILTERBANK: Signal Processing Toolbox (version 6.2.1, or higher) not installed: no filterbank decomposition.']);
wolffd@0:             b = x;
wolffd@0:         end
wolffd@0:     elseif strcmpi(option.filtertype,'Manual');
wolffd@0:         output = cell(1,length(d));
wolffd@0:         for i = 1:length(d)
wolffd@0:             freqi = option.freq;
wolffd@0:             j = 1;
wolffd@0:             while j <= length(freqi)
wolffd@0:                 if not(isinf(freqi(j))) && freqi(j)>f{i}/2
wolffd@0:                     if j == length(freqi)
wolffd@0:                         freqi(j) = Inf;
wolffd@0:                     else
wolffd@0:                         freqi(j) = [];
wolffd@0:                         j = j-1;
wolffd@0:                     end
wolffd@0:                 end
wolffd@0:                 j = j+1;
wolffd@0:             end
wolffd@0:             output{i} = cell(1,length(d{i}));
wolffd@0:             step = option.overlap;
wolffd@0:             for j = 1:length(freqi)-step
wolffd@0:                 if isinf(freqi(j))
wolffd@0:                     [z{j},p{j},k{j}] = ellip(option.filterorder,3,40,...
wolffd@0:                                     freqi(j+step)/f{i}*2);
wolffd@0:                 elseif isinf(freqi(j+step))
wolffd@0:                     [z{j},p{j},k{j}] = ellip(option.filterorder,3,40,...
wolffd@0:                                     freqi(j)/f{i}*2,'high');
wolffd@0:                 else
wolffd@0:                     [z{j},p{j},k{j}] = ellip(option.filterorder,3,40,...
wolffd@0:                                     freqi([j j+step])/f{i}*2);
wolffd@0:                 end
wolffd@0:             end
wolffd@0:             for j = 1:length(z)
wolffd@0:                 [sos,g] = zp2sos(z{j},p{j},k{j});
wolffd@0:                 Hd(j) = dfilt.df2tsos(sos,g);
wolffd@0:                 Hd(j).PersistentMemory = true;
wolffd@0:                 for h = 1:length(d{i})
wolffd@0:                     output{i}{h}(:,:,j) = filter(Hd(j),d{i}{h});
wolffd@0:                 end
wolffd@0:             end
wolffd@0:             %fvtool(Hd)
wolffd@0:             nch{i} = 1:length(freqi)-step;
wolffd@0:         end
wolffd@0:         b = set(x,'Data',output,'Channels',nch);
wolffd@0:     end
wolffd@0:  end
wolffd@0: 
wolffd@0:     
wolffd@0: function [output tmp] = ERBfilterbank(x, fcoefs, chans, tmp)
wolffd@0: % The following is based on the source code from Auditory Toolbox 
wolffd@0: % (A part that I could not call directly from MIRtoolbox, and that I
wolffd@0: % generalize to chunk decomposition)
wolffd@0: 
wolffd@0: % (Malcolm Slaney, August 1993, (c) 1998 Interval Research Corporation)
wolffd@0: 
wolffd@0: % Process an input waveform with a gammatone filter bank.
wolffd@0: % The fcoefs parameter, which completely specifies the Gammatone filterbank,
wolffd@0: % should be designed with the MakeERBFilters function.
wolffd@0: 
wolffd@0: A0  = fcoefs(:,1);
wolffd@0: A11 = fcoefs(:,2);
wolffd@0: A12 = fcoefs(:,3);
wolffd@0: A13 = fcoefs(:,4);
wolffd@0: A14 = fcoefs(:,5);
wolffd@0: A2  = fcoefs(:,6);
wolffd@0: B0  = fcoefs(:,7);
wolffd@0: B1  = fcoefs(:,8);
wolffd@0: B2  = fcoefs(:,9);
wolffd@0: gain= fcoefs(:,10);	
wolffd@0: lc = length(chans);
wolffd@0: output = zeros(size(x,1),size(x,2),lc);
wolffd@0: if isempty(tmp)
wolffd@0:     emptytmp = 1;
wolffd@0: else
wolffd@0:     emptytmp = 0;
wolffd@0: end
wolffd@0: for i = 1:lc
wolffd@0:     chan = length(gain)-chans(i)+1; % Revert the channels order
wolffd@0:     aa1 = [A0(chan)/gain(chan) A11(chan)/gain(chan)  A2(chan)/gain(chan)];
wolffd@0:     bb1 = [B0(chan) B1(chan) B2(chan)];
wolffd@0:     aa2 = [A0(chan) A12(chan) A2(chan)];
wolffd@0:     bb2 = [B0(chan) B1(chan) B2(chan)];
wolffd@0:     aa3 = [A0(chan) A12(chan) A2(chan)];
wolffd@0:     bb3 = [B0(chan) B1(chan) B2(chan)];
wolffd@0:     aa4 = [A0(chan) A13(chan) A2(chan)];
wolffd@0:     bb4 = [B0(chan) B1(chan) B2(chan)];
wolffd@0:     if emptytmp
wolffd@0:         [y1 tmp(:,:,i,1)] = filter(aa1,bb1,x);
wolffd@0:         [y2 tmp(:,:,i,2)] = filter(aa2,bb2,y1);
wolffd@0:         [y3 tmp(:,:,i,3)] = filter(aa3,bb3,y2);
wolffd@0:         [y4 tmp(:,:,i,4)] = filter(aa4,bb4,y3);
wolffd@0:     else
wolffd@0:         [y1 tmp(:,:,i,1)] = filter(aa1,bb1,x,tmp(:,:,i,1));
wolffd@0:         [y2 tmp(:,:,i,2)] = filter(aa2,bb2,y1,tmp(:,:,i,2));
wolffd@0:         [y3 tmp(:,:,i,3)] = filter(aa3,bb3,y2,tmp(:,:,i,3));
wolffd@0:         [y4 tmp(:,:,i,4)] = filter(aa4,bb4,y3,tmp(:,:,i,4));
wolffd@0:     end
wolffd@0:     output(:,:,i) = y4;
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: function [y orig] = eachchunk(orig,option,missing)
wolffd@0: y = mirfilterbank(orig,option);
wolffd@0: 
wolffd@0: 
wolffd@0: function y = combinechunk(old,new)
wolffd@0: do = get(old,'Data');
wolffd@0: to = get(old,'Time');
wolffd@0: dn = get(new,'Data');
wolffd@0: tn = get(new,'Time');
wolffd@0: y = set(old,'Data',{{[do{1}{1};dn{1}{1}]}},...
wolffd@0:             'Time',{{[to{1}{1};tn{1}{1}]}});