wolffd@0: function varargout = mirspectrum(orig,varargin)
wolffd@0: %   s = mirspectrum(x) computes the spectrum of the audio signal x, showing
wolffd@0: %       the distribution of the energy along the frequencies.
wolffd@0: %       (x can be the name of an audio file as well.)
wolffd@0: %   Optional argument:
wolffd@0: %       mirspectrum(...,'Frame',l,h) computes spectrogram, using frames of
wolffd@0: %           length l seconds and a hop rate h.
wolffd@0: %           Default values: l = .05 s, h = .5.
wolffd@0: %       mirspectrum(...,'Min',mi) indicates the lowest frequency taken into
wolffd@0: %           consideration, expressed in Hz.
wolffd@0: %           Default value: 0 Hz.
wolffd@0: %       mirspectrum(...,'Max',ma) indicates the highest frequency taken into
wolffd@0: %           consideration, expressed in Hz.
wolffd@0: %           Default value: the maximal possible frequency, corresponding to
wolffd@0: %           the sampling rate of x divided by 2.
wolffd@0: %       mirspectrum(...,'Window',w): windowing method:
wolffd@0: %           either w = 0 (no windowing) or any windowing function proposed
wolffd@0: %               in the Signal Processing Toolbox (help window).
wolffd@0: %           default value: w = 'hamming'
wolffd@0: %
wolffd@0: %       mirspectrum(...,'Cents'): decomposes the energy in cents.
wolffd@0: %       mirspectrum(...,'Collapsed'): collapses the spectrum into one
wolffd@0: %           octave divided into 1200 cents.
wolffd@0: %       Redistribution of the frequencies into bands:
wolffd@0: %           mirspectrum(...,'Mel'): Mel bands.
wolffd@0: %               (Requires the Auditory Toolbox.)
wolffd@0: %           mirspectrum(...,'Bark'): Bark bands.
wolffd@0: %               (Code based on Pampalk's MA toolbox).
wolffd@0: %           If the audio signal was frame decomposed, the output s is a
wolffd@0: %               band-decomposed spectrogram. It is then possible to compute
wolffd@0: %               the spectrum of the temporal signal in each band,
wolffd@0: %               using the following syntax:
wolffd@0: %                   mirspectrum(s,'AlongBands')
wolffd@0: %               This corresponds to fluctuation (cf. mirfluctuation).
wolffd@0: %           mirspectrum(...,'Mask'): Models masking phenomena in each band.
wolffd@0: %               (Code based on Pampalk's MA toolbox).
wolffd@0: %       mirspectrum(...,'Normal'): normalizes with respect to energy.
wolffd@0: %       mirspectrum(...,'NormalInput'): input signal is normalized from 0 to 1.
wolffd@0: %       mirspectrum(...,'Power'): squares the energy.
wolffd@0: %       mirspectrum(...,'dB'): represents the spectrum energy in decibel scale.
wolffd@0: %       mirspectrum(...,'dB',th): keeps only the highest energy over a
wolffd@0: %           range of th dB.
wolffd@0: %       mirspectrum(...,'Terhardt'): modulates the energy following
wolffd@0: %           (Terhardt, 1979) outer ear model. (Code based on Pampalk's MA
wolffd@0: %           toolbox).
wolffd@0: %       mirspectrum(...,'Resonance',r): multiplies the spectrum with a 
wolffd@0: %           resonance curve. Possible value for r:
wolffd@0: %               'ToiviainenSnyder': highlights best perceived meter 
wolffd@0: %                       (Toiviainen & Snyder 2003)
wolffd@0: %                   (default value for spectrum of envelopes)
wolffd@0: %               'Fluctuation': fluctuation strength (Fastl 1982)
wolffd@0: %                   (default value for spectrum of band channels)
wolffd@0: %       mirspectrum(...,'Prod',m): Enhances components that have harmonics
wolffd@0: %           located at multiples of range(s) m of the signal's fundamental 
wolffd@0: %           frequency. Computed by compressing the signal by the list of
wolffd@0: %           factors m, and by multiplying all the results with the original 
wolffd@0: %           signa (Alonso et al, 2003).
wolffd@0: %       mirspectrum(...,'Sum',s): Similar option using summation instead of
wolffd@0: %           product.
wolffd@0: %
wolffd@0: %       mirspectrum(...,'MinRes',mr): Indicates a minimal accepted
wolffd@0: %           frequency resolution, in Hz. The audio signal is zero-padded in
wolffd@0: %           order to reach the desired resolution.
wolffd@0: %           If the 'Mel' option is toggled on, 'MinRes' is set by default
wolffd@0: %               to 66 Hz.
wolffd@0: %       mirspectrum(...,'MinRes',mr,'OctaveRatio',tol): Indicates the
wolffd@0: %           minimal accepted resolution in terms of number of divisions of 
wolffd@0: %           the octave. Low  frequencies are ignored in order to reach the
wolffd@0: %           desired resolution.
wolffd@0: %               The corresponding required frequency resolution is equal to
wolffd@0: %               the difference between the first frequency bins, multiplied
wolffd@0: %               by the constraining multiplicative factor tol (set by
wolffd@0: %               default to .75).
wolffd@0: %       mirspectrum(...,'Res',r): Indicates the required precise frequency
wolffd@0: %           resolution, in Hz. The audio signal is zero-padded in order to
wolffd@0: %           reach the desired resolution.
wolffd@0: %       mirspectrum(...,'Length',l): Specifies the length of the FFT,
wolffd@0: %           overriding the FFT length initially planned.
wolffd@0: %       mirspectrum(...,'ZeroPad',s): Zero-padding of s samples.
wolffd@0: %       mirspectrum(...,'WarningRes',s): Indicates a required frequency
wolffd@0: %           resolution, in Hz, for the input signal. If the resolution does
wolffd@0: %           not reach that prerequisite, a warning is displayed.
wolffd@0: %       mirspectrum(...,'ConstantQ',nb): Carries out a Constant Q Transform
wolffd@0: %           instead of a FFT, with a number of bins per octave fixed to nb.
wolffd@0: %           Default value for nb: 12 bins per octave.
wolffd@0: %
wolffd@0: %       mirspectrum(...,'Smooth',o): smooths the envelope using a movering
wolffd@0: %           average of order o.
wolffd@0: %           Default value when the option is toggled on: o=10
wolffd@0: %       mirspectrum(...,'Gauss',o): smooths the envelope using a gaussian
wolffd@0: %           of standard deviation o samples.
wolffd@0: %           Default value when the option is toggled on: o=10
wolffd@0: %       mirspectrum(...,'Phase',0): do not compute the FFT phase.
wolffd@0:     
wolffd@0:         win.key = 'Window';
wolffd@0:         win.type = 'String';
wolffd@0:         win.default = 'hamming';
wolffd@0:     option.win = win;
wolffd@0:     
wolffd@0:         min.key = 'Min';
wolffd@0:         min.type = 'Integer';
wolffd@0:         min.default = 0;
wolffd@0:     option.min = min;
wolffd@0:     
wolffd@0:         max.key = 'Max';
wolffd@0:         max.type = 'Integer';
wolffd@0:         max.default = Inf;
wolffd@0:     option.max = max;
wolffd@0:     
wolffd@0:         mr.key = 'MinRes';
wolffd@0:         mr.type = 'Integer';
wolffd@0:         mr.default = 0;
wolffd@0:     option.mr = mr;
wolffd@0: 
wolffd@0:         res.key = 'Res';
wolffd@0:         res.type = 'Integer';
wolffd@0:         res.default = NaN;
wolffd@0:     option.res = res;
wolffd@0: 
wolffd@0:         length.key = 'Length';
wolffd@0:         length.type = 'Integer';
wolffd@0:         length.default = NaN;
wolffd@0:     option.length = length;
wolffd@0:     
wolffd@0:         zp.key = 'ZeroPad';
wolffd@0:         zp.type = 'Integer';
wolffd@0:         zp.default = 0;
wolffd@0:         zp.keydefault = Inf;
wolffd@0:     option.zp = zp;
wolffd@0:     
wolffd@0:         wr.key = 'WarningRes';
wolffd@0:         wr.type = 'Integer';
wolffd@0:         wr.default = 0;
wolffd@0:     option.wr = wr;
wolffd@0:     
wolffd@0:         octave.key = 'OctaveRatio';
wolffd@0:         octave.type = 'Boolean';
wolffd@0:         octave.default = 0;
wolffd@0:     option.octave = octave;
wolffd@0:     
wolffd@0:         constq.key = 'ConstantQ';
wolffd@0:         constq.type = 'Integer';
wolffd@0:         constq.default = 0;
wolffd@0:         constq.keydefault = 12;
wolffd@0:     option.constq = constq;
wolffd@0:     
wolffd@0:         alongbands.key = 'AlongBands';
wolffd@0:         alongbands.type = 'Boolean';
wolffd@0:         alongbands.default = 0;
wolffd@0:     option.alongbands = alongbands;
wolffd@0: 
wolffd@0:         ni.key = 'NormalInput';
wolffd@0:         ni.type = 'Boolean';
wolffd@0:         ni.default = 0;
wolffd@0:     option.ni = ni;
wolffd@0:     
wolffd@0:         norm.key = 'Normal';
wolffd@0:         norm.type = 'Integer';
wolffd@0:         norm.default = 0;
wolffd@0:         norm.keydefault = 1;
wolffd@0:         norm.when = 'After';
wolffd@0:     option.norm = norm;
wolffd@0:     
wolffd@0:         band.type = 'String';
wolffd@0:         band.choice = {'Freq','Mel','Bark','Cents'};
wolffd@0:         band.default = 'Freq';
wolffd@0:         band.when = 'Both';
wolffd@0:     option.band = band;
wolffd@0: 
wolffd@0:         nbbands.key = 'Bands';
wolffd@0:         nbbands.type = 'Integer';
wolffd@0:         nbbands.default = 0;
wolffd@0:         nbbands.when = 'After';
wolffd@0:     option.nbbands = nbbands;
wolffd@0: 
wolffd@0:         mprod.key = 'Prod';
wolffd@0:         mprod.type = 'Integers';
wolffd@0:         mprod.default = [];
wolffd@0:         mprod.keydefault = 2:6;
wolffd@0:         mprod.when = 'After';
wolffd@0:     option.mprod = mprod;
wolffd@0: 
wolffd@0:         msum.key = 'Sum';
wolffd@0:         msum.type = 'Integers';
wolffd@0:         msum.default = [];
wolffd@0:         msum.keydefault = 2:6;
wolffd@0:         msum.when = 'After';
wolffd@0:     option.msum = msum;
wolffd@0: 
wolffd@0:         reso.key = 'Resonance';
wolffd@0:         reso.type = 'String';
wolffd@0:         reso.choice = {'ToiviainenSnyder','Fluctuation','Meter',0,'no','off'};
wolffd@0:         reso.default = 0;
wolffd@0:         reso.when = 'After';
wolffd@0:     option.reso = reso;
wolffd@0:     
wolffd@0:         log.key = 'log';
wolffd@0:         log.type = 'Boolean';
wolffd@0:         log.default = 0;
wolffd@0:         log.when = 'After';
wolffd@0:     option.log = log;
wolffd@0: 
wolffd@0:         db.key = 'dB';
wolffd@0:         db.type = 'Integer';
wolffd@0:         db.default = 0;
wolffd@0:         db.keydefault = Inf;
wolffd@0:         db.when = 'After';
wolffd@0:     option.db = db;
wolffd@0:     
wolffd@0:         pow.key = 'Power';
wolffd@0:         pow.type = 'Boolean';
wolffd@0:         pow.default = 0;
wolffd@0:         pow.when = 'After';
wolffd@0:     option.pow = pow;
wolffd@0: 
wolffd@0:         terhardt.key = 'Terhardt';
wolffd@0:         terhardt.type = 'Boolean';
wolffd@0:         terhardt.default = 0;
wolffd@0:         terhardt.when = 'After';
wolffd@0:     option.terhardt = terhardt;
wolffd@0: 
wolffd@0:         mask.key = 'Mask';
wolffd@0:         mask.type = 'Boolean';
wolffd@0:         mask.default = 0;
wolffd@0:         mask.when = 'After';
wolffd@0:     option.mask = mask;
wolffd@0: 
wolffd@0:     %    e.key = 'Enhanced';
wolffd@0:     %    e.type = 'Integer';
wolffd@0:     %    e.default = [];
wolffd@0:     %    e.keydefault = 2:10;
wolffd@0:     %    e.when = 'After';
wolffd@0:     %option.e = e;
wolffd@0: 
wolffd@0:         collapsed.key = 'Collapsed';
wolffd@0:         collapsed.type = 'Boolean';
wolffd@0:         collapsed.default = 0;
wolffd@0:         collapsed.when = 'Both';
wolffd@0:     option.collapsed = collapsed;
wolffd@0:     
wolffd@0:         aver.key = 'Smooth';
wolffd@0:         aver.type = 'Integer';
wolffd@0:         aver.default = 0;
wolffd@0:         aver.keydefault = 10;
wolffd@0:         aver.when = 'After';
wolffd@0:     option.aver = aver;
wolffd@0:     
wolffd@0:         gauss.key = 'Gauss';
wolffd@0:         gauss.type = 'Integer';
wolffd@0:         gauss.default = 0;
wolffd@0:         gauss.keydefault = 10;
wolffd@0:         gauss.when = 'After';
wolffd@0:     option.gauss = gauss;
wolffd@0:     
wolffd@0:         rapid.key = 'Rapid';
wolffd@0:         rapid.type = 'Boolean';
wolffd@0:         rapid.default = 0;
wolffd@0:     option.rapid = rapid;
wolffd@0: 
wolffd@0:         phase.key = 'Phase';
wolffd@0:         phase.type = 'Boolean';
wolffd@0:         phase.default = 1;
wolffd@0:     option.phase = phase;
wolffd@0: 
wolffd@0: specif.option = option;
wolffd@0: 
wolffd@0: specif.defaultframelength = 0.05;
wolffd@0: specif.defaultframehop = 0.5;
wolffd@0: 
wolffd@0: specif.eachchunk = @eachchunk;
wolffd@0: specif.combinechunk = @combinechunk;
wolffd@0: 
wolffd@0: if isamir(orig,'mirscalar') || isamir(orig,'mirenvelope')
wolffd@0:     specif.nochunk = 1;
wolffd@0: end
wolffd@0: 
wolffd@0: varargout = mirfunction(@mirspectrum,orig,varargin,nargout,specif,@init,@main);
wolffd@0: 
wolffd@0: 
wolffd@0: function [x type] = init(x,option)
wolffd@0: type = 'mirspectrum';
wolffd@0: 
wolffd@0: 
wolffd@0: function s = main(orig,option,postoption)
wolffd@0: if isstruct(option)
wolffd@0:     if option.collapsed
wolffd@0:         option.band = 'Cents';
wolffd@0:     end
wolffd@0:     if isnan(option.res) && strcmpi(option.band,'Cents') && option.min
wolffd@0:         option.res = option.min *(2^(1/1200)-1)*.9;
wolffd@0:     end    
wolffd@0: end
wolffd@0: if not(isempty(postoption))
wolffd@0:     if not(strcmpi(postoption.band,'Freq') && isempty(postoption.msum) ...
wolffd@0:             || isempty(postoption.mprod)) || postoption.log || postoption.db ...
wolffd@0:             || postoption.pow || postoption.mask || postoption.collapsed ...
wolffd@0:             || postoption.aver || postoption.gauss
wolffd@0:         option.phase = 0;
wolffd@0:     end
wolffd@0: end
wolffd@0: if iscell(orig)
wolffd@0:     orig = orig{1};
wolffd@0: end
wolffd@0: if isa(orig,'mirspectrum') && ...
wolffd@0:         not(isfield(option,'alongbands') && option.alongbands)
wolffd@0:     s = orig;
wolffd@0:     if isfield(option,'min') && ...
wolffd@0:             (option.min || iscell(option.max) || option.max < Inf)
wolffd@0:         magn = get(s,'Magnitude');
wolffd@0:         freq = get(s,'Frequency');
wolffd@0:         for k = 1:length(magn)
wolffd@0:             m = magn{k};
wolffd@0:             f = freq{k};
wolffd@0:             if iscell(option.max)
wolffd@0:                 mi = option.min{k};
wolffd@0:                 ma = option.max{k};
wolffd@0:             else
wolffd@0:                 mi = option.min;
wolffd@0:                 ma = option.max;
wolffd@0:             end
wolffd@0:             if not(iscell(m))
wolffd@0:                 m = {m};
wolffd@0:                 f = {f};
wolffd@0:             end
wolffd@0:             for l = 1:length(m)
wolffd@0:                 range = find(and(f{l}(:,1) >= mi,f{l}(:,1) <= ma));
wolffd@0:                 m{l} = m{l}(range,:,:);
wolffd@0:                 f{l} = f{l}(range,:,:);
wolffd@0:             end
wolffd@0:             magn{k} = m;
wolffd@0:             freq{k} = f;
wolffd@0:         end
wolffd@0:         s = set(s,'Magnitude',magn,'Frequency',freq);
wolffd@0:     end
wolffd@0:     if not(isempty(postoption)) && not(isequal(postoption,0))
wolffd@0:         s = post(s,postoption);
wolffd@0:     end
wolffd@0: elseif ischar(orig)
wolffd@0:     s = mirspectrum(miraudio(orig),option,postoption);
wolffd@0: else
wolffd@0:     if nargin == 0
wolffd@0:         orig = [];
wolffd@0:     end
wolffd@0:     s.phase = [];
wolffd@0:     s.log = 0;
wolffd@0:     s.xscale = 'Freq';
wolffd@0:     s.pow = 1;
wolffd@0:     s = class(s,'mirspectrum',mirdata(orig));
wolffd@0:     s = purgedata(s);
wolffd@0:     s = set(s,'Title','Spectrum','Abs','frequency (Hz)','Ord','magnitude');
wolffd@0:     %disp('Computing spectrum...')
wolffd@0:     sig = get(orig,'Data');
wolffd@0:     t = get(orig,'Pos');
wolffd@0:     if isempty(t)
wolffd@0:         t = get(orig,'FramePos');
wolffd@0:         for k = 1:length(sig)
wolffd@0:             for l = 1:length(sig{k})
wolffd@0:                 sig{k}{l} = sig{k}{l}';
wolffd@0:                 t{k}{l} = t{k}{l}(1,:,:)';
wolffd@0:             end
wolffd@0:         end
wolffd@0:     end
wolffd@0:     fs = get(orig,'Sampling');
wolffd@0:     fp = get(orig,'FramePos');
wolffd@0:     m = cell(1,length(sig));
wolffd@0:     p = cell(1,length(sig));
wolffd@0:     f = cell(1,length(sig));
wolffd@0:     for i = 1:length(sig)
wolffd@0:         d = sig{i};
wolffd@0:         fpi = fp{i};
wolffd@0:         if not(iscell(d))
wolffd@0:             d = {d};
wolffd@0:         end
wolffd@0:         if option.alongbands
wolffd@0:             fsi = 1 / (fpi{1}(1,2) - fpi{1}(1,1));
wolffd@0:         else
wolffd@0:             fsi = fs{i};
wolffd@0:         end
wolffd@0:         mi = cell(1,length(d));
wolffd@0:         phi = cell(1,length(d));
wolffd@0:         fi = cell(1,length(d));
wolffd@0:         for J = 1:length(d)
wolffd@0:             dj = d{J};
wolffd@0:             if option.ni
wolffd@0:                 mxdj = repmat(max(dj),[size(dj,1),1,1]);
wolffd@0:                 mndj = repmat(min(dj),[size(dj,1),1,1]);
wolffd@0:                 dj = (dj-mndj)./(mxdj-mndj);
wolffd@0:             end
wolffd@0:             if option.alongbands
wolffd@0:                 if size(dj,1)>1
wolffd@0:                     error('ERROR IN MIRSPECTRUM: ''AlongBands'' is restricted to spectrum decomposed into bands, such as ''Mel'' and ''Bark''.') 
wolffd@0:                 end
wolffd@0:                 dj = reshape(dj,[size(dj,2),1,size(dj,3)]);
wolffd@0:                 fp{i}{J} = fp{i}{J}([1;end]);
wolffd@0:             end
wolffd@0:                         
wolffd@0:             if option.constq
wolffd@0:                 % Constant Q Transform
wolffd@0:                 r = 2^(1/option.constq);
wolffd@0:                 Q = 1 / (r - 1);
wolffd@0:                 f_max = min(fsi/2,option.max);
wolffd@0:                 f_min = option.min;
wolffd@0:                 if not(f_min)
wolffd@0:                     f_min = 16.3516;
wolffd@0:                 end
wolffd@0:                 B = floor(log(f_max/f_min) / log(r)); % number of bins
wolffd@0:                 N0 = round(Q*fsi/f_min); % maximum Nkcq
wolffd@0:                 j2piQn = -j*2*pi*Q*(0:N0-1)';
wolffd@0: 
wolffd@0:                 fj = f_min * r.^(0:B-1)';
wolffd@0:                 transf = NaN(B,size(dj,2),size(dj,3));
wolffd@0:                 for kcq = 1:B
wolffd@0:                     Nkcq = round(Q*fsi/fj(kcq));
wolffd@0:                     win = mirwindow(dj,option.win,Nkcq);
wolffd@0:                     exq = repmat(exp(j2piQn(1:Nkcq)/Nkcq),...
wolffd@0:                                  [1,size(win,2),size(win,3)]);
wolffd@0:                     transf(kcq,:) = sum(win.* exq) / Nkcq;
wolffd@0:                 end
wolffd@0:             else
wolffd@0:                 % FFT
wolffd@0:                 dj = mirwindow(dj,option.win);
wolffd@0:                 if option.zp
wolffd@0:                     if option.zp < Inf
wolffd@0:                         dj = [dj;zeros(option.zp,size(dj,2),size(dj,3))];
wolffd@0:                     else
wolffd@0:                         dj = [dj;zeros(size(dj))];
wolffd@0:                     end
wolffd@0:                 end
wolffd@0:                 if isstruct(postoption)
wolffd@0:                     if strcmpi(postoption.band,'Mel') && ...
wolffd@0:                             (not(option.mr) || option.mr > 66)
wolffd@0:                         option.mr = 66;
wolffd@0:                     end
wolffd@0:                 else
wolffd@0:                     %warning('WARNING in MIRSPECTRUM (for debugging purposes): By default, minimum resolution specified.')
wolffd@0:                     if not(option.mr)
wolffd@0:                         option.mr = 66;
wolffd@0:                     end
wolffd@0:                 end
wolffd@0:                 if option.octave
wolffd@0:                     N = size(dj,1);
wolffd@0:                     res = (2.^(1/option.mr)-1)*option.octave;
wolffd@0:                         % Minimal freq ratio between 2 first bins.
wolffd@0:                         % freq resolution should be > option.min * res
wolffd@0:                     Nrec = fsi/(option.min*res);
wolffd@0:                         % Recommended minimal sample length.
wolffd@0:                     if Nrec > N
wolffd@0:                             % If actual sample length is too small.
wolffd@0:                         option.min = fsi/N / res;
wolffd@0:                         warning('WARNING IN MIRSPECTRUM: The input signal is too short to obtain the desired octave resolution. Lowest frequencies will be ignored.');
wolffd@0:                         display(['(The recommended minimal input signal length would be ' num2str(Nrec/fsi) ' s.)']);
wolffd@0:                         display(['New low frequency range: ' num2str(option.min) ' Hz.']);
wolffd@0:                     end
wolffd@0:                     N = 2^nextpow2(N);
wolffd@0:                 elseif isnan(option.length)
wolffd@0:                     if isnan(option.res)
wolffd@0:                         N = size(dj,1);
wolffd@0:                         if option.mr && N < fsi/option.mr
wolffd@0:                             if option.wr && N < fsi/option.wr
wolffd@0:                                 warning('WARNING IN MIRSPECTRUM: The input signal is too short to obtain the desired frequency resolution. Performed zero-padding will not guarantee the quality of the results.'); 
wolffd@0:                             end                
wolffd@0:                             N = max(N,fsi/option.mr);
wolffd@0:                         end
wolffd@0:                         N = 2^nextpow2(N);
wolffd@0:                     else
wolffd@0:                         N = ceil(fsi/option.res);
wolffd@0:                     end
wolffd@0:                 else
wolffd@0:                     N = option.length;
wolffd@0:                 end
wolffd@0: 
wolffd@0:                 % Here is the spectrum computation itself
wolffd@0:                 transf = fft(dj,N);
wolffd@0: 
wolffd@0:                 len = ceil(size(transf,1)/2);
wolffd@0:                 fj = fsi/2 * linspace(0,1,len)';
wolffd@0:                 if option.max
wolffd@0:                     maxf = find(fj>=option.max,1);
wolffd@0:                     if isempty(maxf)
wolffd@0:                         maxf = len;
wolffd@0:                     end
wolffd@0:                 else
wolffd@0:                     maxf = len;
wolffd@0:                 end
wolffd@0:                 if option.min
wolffd@0:                     minf = find(fj>=option.min,1);
wolffd@0:                     if isempty(minf)
wolffd@0:                         maxf = len;
wolffd@0:                     end
wolffd@0:                 else
wolffd@0:                     minf = 1;
wolffd@0:                 end
wolffd@0:                 
wolffd@0:                 transf = transf(minf:maxf,:,:);
wolffd@0:                 fj = fj(minf:maxf);
wolffd@0:             end
wolffd@0:             
wolffd@0:             mi{J} = abs(transf);
wolffd@0:             if option.phase
wolffd@0:                 phi{J} = angle(transf);
wolffd@0:             end
wolffd@0:             fi{J} = repmat(fj,[1,size(transf,2),size(transf,3)]);
wolffd@0:         end
wolffd@0:         if iscell(sig{i})
wolffd@0:             m{i} = mi;
wolffd@0:             p{i} = phi;
wolffd@0:             f{i} = fi;
wolffd@0:         else
wolffd@0:             m{i} = mi{1};
wolffd@0:             p{i} = phi{1};
wolffd@0:             f{i} = fi{1};
wolffd@0:         end
wolffd@0:     end
wolffd@0:     s = set(s,'Frequency',f,'Magnitude',m,'Phase',p,'FramePos',fp);
wolffd@0:     if not(isempty(postoption)) && isstruct(postoption)
wolffd@0:         s = post(s,postoption);
wolffd@0:     end
wolffd@0: end
wolffd@0: 
wolffd@0:    
wolffd@0: function s = post(s,option)
wolffd@0: if option.collapsed
wolffd@0:     option.band = 'Cents';
wolffd@0: end
wolffd@0: m = get(s,'Magnitude');
wolffd@0: f = get(s,'Frequency');
wolffd@0: for k = 1:length(m)
wolffd@0:     if not(iscell(m{k}))
wolffd@0:         m{k} = {m{k}};
wolffd@0:         f{k} = {f{k}};
wolffd@0:     end
wolffd@0: end
wolffd@0: if get(s,'Power') == 1 && ...
wolffd@0:         (option.pow || any(option.mprod) || any(option.msum)) 
wolffd@0:     for h = 1:length(m)
wolffd@0:         for l = 1:length(m{k})
wolffd@0:             m{h}{l} = m{h}{l}.^2;
wolffd@0:         end
wolffd@0:     end
wolffd@0:     s = set(s,'Power',2,'Title',['Power ',get(s,'Title')]);
wolffd@0: end
wolffd@0: if any(option.mprod)
wolffd@0:     for h = 1:length(m)
wolffd@0:         for l = 1:length(m{k})
wolffd@0:             z0 = m{h}{l};
wolffd@0:             z1 = z0;
wolffd@0:             for k = 1:length(option.mprod)
wolffd@0:                 mpr = option.mprod(k);
wolffd@0:                 if mpr
wolffd@0:                     zi = ones(size(z0));
wolffd@0:                     zi(1:floor(end/mpr),:,:) = z0(mpr:mpr:end,:,:);
wolffd@0:                     z1 = z1.*zi;
wolffd@0:                 end
wolffd@0:             end
wolffd@0:             m{h}{l} = z1;
wolffd@0:         end
wolffd@0:     end
wolffd@0:     s = set(s,'Title','Spectral product');
wolffd@0: end
wolffd@0: if any(option.msum)
wolffd@0:     for h = 1:length(m)
wolffd@0:         for l = 1:length(m{k})
wolffd@0:             z0 = m{h}{l};
wolffd@0:             z1 = z0;
wolffd@0:             for k = 1:length(option.msum)
wolffd@0:                 mpr = option.msum(k);
wolffd@0:                 if mpr
wolffd@0:                     zi = ones(size(z0));
wolffd@0:                     zi(1:floor(end/mpr),:,:) = z0(mpr:mpr:end,:,:);
wolffd@0:                     z1 = z1+zi;
wolffd@0:                 end
wolffd@0:             end
wolffd@0:             m{h}{l} = z1;
wolffd@0:         end
wolffd@0:     end
wolffd@0:     s = set(s,'Title','Spectral sum');
wolffd@0: end
wolffd@0: %for i = option.e  % Enhancing procedure %%%% TO CHECK, t IS NOT DEFINED!
wolffd@0: %    for h = 1:length(m)
wolffd@0: %        for l = 1:length(m{k})
wolffd@0: %            c = m{h}{l};
wolffd@0: %            % option.e is the list of scaling factors
wolffd@0: %            % i is the scaling factor
wolffd@0: %            if i
wolffd@0: %                ic = zeros(size(c));
wolffd@0: %                for g = 1:size(c,2)
wolffd@0: %                    for h2 = 1:size(c,3)
wolffd@0: %                        beg = find(t(:,g,h2)/i >= t(1,g,h2))
wolffd@0: %                        if not(isempty(beg))
wolffd@0: %                            ic(:,g,h2) = interp1(t(:,g,h2),c(:,g,h2),t(:,g,h2)/i); %,'cubic'); 
wolffd@0:                                 % The scaled autocorrelation
wolffd@0: %                            ic(1:beg(1)-1,g,h2) = 0;
wolffd@0: %                        end
wolffd@0: %                    end
wolffd@0: %                end
wolffd@0: %                ic(find(isnan(ic))) = Inf;    % All the NaN values are changed into 0 in the resulting curve
wolffd@0: %                ic = max(ic,0);
wolffd@0: %                c = c - ic;       % The scaled autocorrelation is substracted to the initial one
wolffd@0: %                c = max(c,0);               % Half-wave rectification
wolffd@0:                 %figure
wolffd@0:                 %plot(c)
wolffd@0: %            end
wolffd@0: %        end
wolffd@0: %    end
wolffd@0: %end
wolffd@0: if option.norm
wolffd@0:     for k = 1:length(m)
wolffd@0:         for l = 1:length(m{k})
wolffd@0:             mkl = m{k}{l};
wolffd@0:             nkl = zeros(1,size(mkl,2),size(mkl,3));
wolffd@0:             for kk = 1:size(mkl,2)
wolffd@0:                 for ll = 1:size(mkl,3)
wolffd@0:                     nkl(1,kk,l) = norm(mkl(:,kk,ll));
wolffd@0:                 end
wolffd@0:             end
wolffd@0:             m{k}{l} = mkl./repmat(nkl,[size(m{k}{k},1),1,1]);
wolffd@0:         end
wolffd@0:     end
wolffd@0: end
wolffd@0: if option.terhardt && not(isempty(find(f{1}{1}))) % This excludes the case where spectrum already along bands
wolffd@0:     % Code taken from Pampalk's MA Toolbox
wolffd@0:     for h = 1:length(m)
wolffd@0:         for l = 1:length(m{k})
wolffd@0:             W_Adb = zeros(size(f{k}{l}));
wolffd@0:             W_Adb(2:size(f{k}{l},1),:,:) = ...
wolffd@0:                 + 10.^((-3.64*(f{k}{l}(2:end,:,:)/1000).^-0.8 ...
wolffd@0:                 + 6.5 * exp(-0.6 * (f{k}{l}(2:end,:,:)/1000 - 3.3).^2) ...
wolffd@0:                 - 0.001*(f{k}{l}(2:end,:,:)/1000).^4)/20);
wolffd@0:             W_Adb = W_Adb.^2;
wolffd@0:             m{h}{l} = m{h}{l}.*W_Adb;
wolffd@0:         end
wolffd@0:     end
wolffd@0: end
wolffd@0: if option.reso
wolffd@0:     if not(ischar(option.reso))
wolffd@0:         if strcmp(get(orig,'XScale'),'Mel')
wolffd@0:             option.reso = 'Fluctuation';
wolffd@0:         else
wolffd@0:             option.reso = 'ToiviainenSnyder';
wolffd@0:         end
wolffd@0:     end
wolffd@0:     %if strcmpi(option.reso,'ToiviainenSnyder') || ...
wolffd@0:     %   strcmpi(option.reso,'Meter') || ...
wolffd@0:     %   strcmpi(option.reso,'Fluctuation')
wolffd@0:     for k = 1:length(m)
wolffd@0:         for l = 1:length(m{k})
wolffd@0:             if strcmpi(option.reso,'ToiviainenSnyder') || strcmpi(option.reso,'Meter')
wolffd@0:                 w = max(0,...
wolffd@0:                     1 - 0.25*(log2(max(1./max(f{k}{l},1e-12),1e-12)/0.5)).^2);
wolffd@0:             elseif strcmpi(option.reso,'Fluctuation')
wolffd@0:                 w1 = f{k}{l} / 4; % ascending part of the fluctuation curve;
wolffd@0:                 w2 = 1 - 0.3 * (f{k}{l} - 4)/6; % descending part;
wolffd@0:                 w = min(w1,w2);
wolffd@0:             end
wolffd@0:             if max(w) == 0
wolffd@0:                 warning('The resonance curve, not defined for this range of delays, will not be applied.')
wolffd@0:             else
wolffd@0:                 m{k}{l} = m{k}{l}.* w;
wolffd@0:             end
wolffd@0:         end
wolffd@0:     end
wolffd@0: end
wolffd@0: if strcmp(s.xscale,'Freq')
wolffd@0:     if strcmpi(option.band,'Mel') 
wolffd@0:         % The following is largely based on the source code from Auditory Toolbox 
wolffd@0:         % (A part that I could not call directly from MIRtoolbox)
wolffd@0: 
wolffd@0:         % (Malcolm Slaney, August 1993, (c) 1998 Interval Research Corporation)
wolffd@0: 
wolffd@0:         try
wolffd@0:             MakeERBFilters(1,1,1); % Just to be sure that the Auditory Toolbox is installed
wolffd@0:         catch
wolffd@0:             error('ERROR IN MIRFILTERBANK: Auditory Toolbox needs to be installed.');
wolffd@0:         end  
wolffd@0: 
wolffd@0:         %disp('Computing Mel-frequency spectral representation...')
wolffd@0:         lowestFrequency = 133.3333;
wolffd@0:         if not(option.nbbands)
wolffd@0:             option.nbbands = 40;
wolffd@0:         end
wolffd@0:         linearFilters = min(13,option.nbbands);
wolffd@0:         linearSpacing = 66.66666666;
wolffd@0:         logFilters = option.nbbands - linearFilters;
wolffd@0:         logSpacing = 1.0711703;
wolffd@0:         totalFilters = option.nbbands;
wolffd@0: 
wolffd@0:         % Figure the band edges.  Interesting frequencies are spaced
wolffd@0:         % by linearSpacing for a while, then go logarithmic.  First figure
wolffd@0:         % all the interesting frequencies.  Lower, center, and upper band
wolffd@0:         % edges are all consequtive interesting frequencies. 
wolffd@0:         freqs = lowestFrequency + (0:linearFilters-1)*linearSpacing;
wolffd@0:         freqs(linearFilters+1:totalFilters+2) = ...
wolffd@0:                       freqs(linearFilters) * logSpacing.^(1:logFilters+2);
wolffd@0:         lower = freqs(1:totalFilters);
wolffd@0:         center = freqs(2:totalFilters+1);
wolffd@0:         upper = freqs(3:totalFilters+2);
wolffd@0:         
wolffd@0:         % Figure out the height of the triangle so that each filter has 
wolffd@0:         % unit weight, assuming a triangular weighting function.
wolffd@0:         triangleHeight = 2./(upper-lower);
wolffd@0: 
wolffd@0:         e = cell(1,length(m));
wolffd@0:         nch = cell(1,length(m));
wolffd@0:         for h = 1:length(m)
wolffd@0:             e{h} = cell(1,length(m{h}));
wolffd@0:             for i = 1:length(m{h})
wolffd@0:                 mi = m{h}{i};
wolffd@0:                 fi = f{h}{i}(:,1,1);
wolffd@0:                 fftSize = size(fi,1);
wolffd@0:                 
wolffd@0:                 % We now want to combine FFT bins and figure out 
wolffd@0:                 % each frequencies contribution
wolffd@0:                 mfccFilterWeights = zeros(totalFilters,fftSize);
wolffd@0:                 for chan=1:totalFilters
wolffd@0:                      mfccFilterWeights(chan,:) = triangleHeight(chan).*...
wolffd@0:                         ((fi > lower(chan) & fi <= center(chan)).* ...
wolffd@0:                          (fi-lower(chan))/(center(chan)-lower(chan)) + ...
wolffd@0:                          (fi > center(chan) & fi < upper(chan)).* ...
wolffd@0:                          (upper(chan)-fi)/(upper(chan)-center(chan)));
wolffd@0:                 end
wolffd@0:                 if find(diff(not(sum(mfccFilterWeights,2)))==-1)
wolffd@0:                     % If one channel has no weight whereas the higher one
wolffd@0:                     % has a positive weight.
wolffd@0:                     warning('WARNING in MIRSPECTRUM: The frequency resolution of the spectrum is too low for the Mel transformation. Some Mel components are undefined.')
wolffd@0:                     display('Recommended frequency resolution: at least 66 Hz.')
wolffd@0:                 end
wolffd@0:                 e{h}{i} = zeros(1,size(mi,2),totalFilters);
wolffd@0:                 for J = 1:size(mi,2)
wolffd@0:                     if max(mi(:,J)) > 0
wolffd@0:                         fftData = zeros(fftSize,1); % Zero-padding ?
wolffd@0:                         fftData(1:size(mi,1)) = mi(:,J);
wolffd@0:                         p = mfccFilterWeights * fftData + 1e-16;
wolffd@0:                         e{h}{i}(1,J,:) = reshape(p,[1,1,length(p)]);
wolffd@0:                     end
wolffd@0:                 end
wolffd@0:                 f{h}{i} = zeros(1,size(mi,2),totalFilters);
wolffd@0:             end
wolffd@0:             nch{h} = (1:totalFilters)';
wolffd@0:         end
wolffd@0:         m = e;
wolffd@0:         s = set(s,'XScale','Mel','Title','Mel-Spectrum','Abs','Mel bands','Channels',nch);
wolffd@0:     elseif strcmpi(option.band,'Bark')
wolffd@0:         sr = get(s,'Sampling');
wolffd@0:         % Code taken from Pampalk's MA Toolbox.
wolffd@0:         %% zwicker & fastl: psychoacoustics 1999, page 159
wolffd@0:         bark_upper = [10 20 30 40 51 63 77 92 108 127 148 172 200 232 270 315 370 440 530 640 770 950 1200 1550]*10; %% Hz
wolffd@0:         e = cell(1,length(m));
wolffd@0:         nch = cell(1,length(m));
wolffd@0:         for h = 1:length(m)
wolffd@0:             %% ignore critical bands outside of sampling rate range
wolffd@0:             cb = min(min([find(bark_upper>sr{h}/2),length(bark_upper)]),length(bark_upper));
wolffd@0:             e{h} = cell(1,length(m{h}));
wolffd@0:             for i = 1:length(m{h})
wolffd@0:                 mi = sum(m{h}{i},3);
wolffd@0:                 e{h}{i} = zeros(1,size(mi,2),cb);
wolffd@0:                 k = 1;
wolffd@0:                 for J=1:cb, %% group into bark bands
wolffd@0:                     idx = find(f{h}{i}(k:end,1,1)<=bark_upper(J));
wolffd@0:                     idx = idx + k-1;
wolffd@0:                     e{h}{i}(1,:,J) = sum(mi(idx,:,:),1);
wolffd@0:                     k = max(idx)+1;
wolffd@0:                 end
wolffd@0:                 f{h}{i} = zeros(1,size(mi,2),cb);
wolffd@0:             end
wolffd@0:             nch{h} = (1:length(bark_upper))';
wolffd@0:            end
wolffd@0:         m = e;
wolffd@0:         s = set(s,'XScale','Bark','Title','Bark-Spectrum','Abs','Bark bands','Channels',nch);
wolffd@0:     elseif strcmpi(option.band,'Cents') || option.collapsed
wolffd@0:         for h = 1:length(m)
wolffd@0:             for i = 1:length(m{h})
wolffd@0:                 mi = m{h}{i};
wolffd@0:                 fi = f{h}{i};
wolffd@0:                 isgood = fi(:,1,1)*(2^(1/1200)-1) >= fi(2,1,1)-fi(1,1,1);
wolffd@0:                 good = find(isgood);
wolffd@0:                 if isempty(good)
wolffd@0:                     mirerror('mirspectrum',...
wolffd@0:                         'The frequency resolution of the spectrum is too low to be decomposed into cents.');
wolffd@0:                     display('Hint: if you specify a minimum value for the frequency range (''Min'' option)');
wolffd@0:                     display('      and if you do not specify any frequency resolution (''Res'' option), ');
wolffd@0:                     display('      the correct frequency resolution will be automatically chosen.');
wolffd@0:                 end
wolffd@0:                 if good>1
wolffd@0:                     warning(['MIRSPECTRUM: Frequency resolution in cent is achieved for frequencies over ',...
wolffd@0:                         num2str(floor(fi(good(1)))),' Hz. Lower frequencies are ignored.'])
wolffd@0:                     display('Hint: if you specify a minimum value for the frequency range (''Min'' option)');
wolffd@0:                     display('      and if you do not specify any frequency resolution (''Res'' option), ');
wolffd@0:                     display('      the correct frequency resolution will be automatically chosen.');
wolffd@0:                 end
wolffd@0:                 fi = fi(good,:,:);
wolffd@0:                 mi = mi(good,:,:);
wolffd@0:                 f2cents = 440*2.^(1/1200*(-1200*10:1200*10-1)');
wolffd@0:                     % The frequencies corresponding to the cent range
wolffd@0:                 cents = repmat((0:1199)',[20,size(fi,2),size(fi,3)]);
wolffd@0:                     % The cent range is for the moment centered to A440
wolffd@0:                 octaves = ones(1200,1)*(-10:10);
wolffd@0:                 octaves = repmat(octaves(:),[1,size(fi,2),size(fi,3)]);
wolffd@0:                 select = find(f2cents>fi(1) & f2cents<fi(end));
wolffd@0:                     % Cent range actually used in the current spectrum
wolffd@0:                 f2cents = repmat(f2cents(select),[1,size(fi,2),size(fi,3)]);
wolffd@0:                 cents = repmat(cents(select),[1,size(fi,2),size(fi,3)]);
wolffd@0:                 octaves = repmat(octaves(select),[1,size(fi,2),size(fi,3)]);
wolffd@0:                 m{h}{i} = interp1(fi(:,1,1),mi,f2cents(:,1,1));
wolffd@0:                 f{h}{i} = octaves*1200+cents + 6900;
wolffd@0:                     % Now the cent range is expressed in midicent
wolffd@0:             end
wolffd@0:         end
wolffd@0:         s = set(s,'Abs','pitch (in midicents)','XScale','Cents');
wolffd@0:     end
wolffd@0: end
wolffd@0: if option.mask
wolffd@0:     if strcmp(s.xscale,'Freq')
wolffd@0:         warning('WARNING IN MIRSPECTRUM: ''Mask'' option available only for Mel-spectrum and Bark-spectrum.');
wolffd@0:         disp('''Mask'' option ignored');
wolffd@0:     else
wolffd@0:         nch = get(s,'Channels');
wolffd@0:         for h = 1:length(m)
wolffd@0:             % Code taken from Pampalk's MA Toolbox.
wolffd@0:             %% spreading function: schroeder et al., 1979, JASA, optimizing
wolffd@0:             %% digital speech coders by exploiting masking properties of the human ear
wolffd@0:             cb = length(nch{h});  % Number of bands.
wolffd@0:             for i = 1:cb, 
wolffd@0:                 spread(i,:) = 10.^((15.81+7.5*((i-(1:cb))+0.474)-17.5*(1+((i-(1:cb))+0.474).^2).^0.5)/10);
wolffd@0:             end
wolffd@0:             for i = 1:length(m{h})
wolffd@0:                 for J = 1:size(m{h}{i},2)
wolffd@0:                     mj = m{h}{i}(1,J,:);
wolffd@0:                     mj = spread(1:length(mj),1:length(mj))*mj(:);
wolffd@0:                     m{h}{i}(1,J,:) = reshape(mj,1,1,length(mj));
wolffd@0:                 end
wolffd@0:             end
wolffd@0:         end
wolffd@0:     end
wolffd@0: end
wolffd@0: if option.collapsed
wolffd@0:     for h = 1:length(m)
wolffd@0:         for i = 1:length(m{h})
wolffd@0:             mi = m{h}{i};
wolffd@0:             fi = f{h}{i};
wolffd@0:             centclass = rem(fi(:,1,1),1200);
wolffd@0:             %neg = find(centclass<0);
wolffd@0:             %centclass(neg) = 1200 + centclass(neg);
wolffd@0:             m{h}{i} = NaN(1200,size(mi,2),size(mi,3));
wolffd@0:             for k = 0:1199
wolffd@0:                 m{h}{i}(k+1,:,:) = sum(mi(find(centclass==k),:,:),1);
wolffd@0:             end
wolffd@0:             f{h}{i} = repmat((0:1199)',[1,size(fi,2),size(fi,3)]);
wolffd@0:         end
wolffd@0:     end
wolffd@0:     s = set(s,'Abs','Cents','XScale','Cents(Collapsed)');
wolffd@0: end
wolffd@0: if option.log || option.db
wolffd@0:     if not(isa(s,'mirspectrum') && s.log)
wolffd@0:         if option.db
wolffd@0:             s = set(s,'log',10);
wolffd@0:         else
wolffd@0:             s = set(s,'log',1);
wolffd@0:         end
wolffd@0:         for k = 1:length(m)
wolffd@0:             if not(iscell(m{k}))
wolffd@0:                 m{k} = {m{k}};
wolffd@0:                 f{k} = {f{k}};
wolffd@0:             end
wolffd@0:             for l = 1:length(m{k})
wolffd@0:                 m{k}{l} = log10(m{k}{l}+1e-16);
wolffd@0:                 if option.db
wolffd@0:                     m{k}{l} = 10*m{k}{l};
wolffd@0:                 end
wolffd@0:             end
wolffd@0:         end
wolffd@0:     elseif isa(s,'mirspectrum') && s.log<10
wolffd@0:         for k = 1:length(m)
wolffd@0:             for l = 1:length(m{k})
wolffd@0:                 m{k}{l} = 10*m{k}{l};
wolffd@0:             end
wolffd@0:         end
wolffd@0:     end
wolffd@0:     if option.db>0 && option.db < Inf
wolffd@0:         for k = 1:length(m)
wolffd@0:             for l = 1:length(m{k})
wolffd@0:                 m{k}{l} = m{k}{l}-repmat(max(m{k}{l}),[size(m{k}{l},1) 1 1]);
wolffd@0:                 m{k}{l} = option.db + max(-option.db,m{k}{l});
wolffd@0:             end
wolffd@0:         end
wolffd@0:     end
wolffd@0: end
wolffd@0: if option.aver
wolffd@0:     for k = 1:length(m)
wolffd@0:         for i = 1:length(m{k})
wolffd@0:             m{k}{i} = filter(ones(1,option.aver),1,m{k}{i});
wolffd@0:         end
wolffd@0:     end
wolffd@0: end
wolffd@0: if option.gauss
wolffd@0:     for k = 1:length(m)
wolffd@0:         for i = 1:length(m{k})
wolffd@0:             sigma = option.gauss;
wolffd@0:             gauss = 1/sigma/2/pi...
wolffd@0:                     *exp(- (-4*sigma:4*sigma).^2 /2/sigma^2);
wolffd@0:             y = filter(gauss,1,[m{k}{i};zeros(4*sigma,1)]);
wolffd@0:             y = y(4*sigma:end,:,:);
wolffd@0:             m{k}{i} = y(1:size(m{k}{i},1),:,:);
wolffd@0:         end
wolffd@0:     end
wolffd@0: end
wolffd@0: s = set(s,'Magnitude',m,'Frequency',f);
wolffd@0: 
wolffd@0: 
wolffd@0: function dj = mirwindow(dj,win,N)
wolffd@0: if nargin<3
wolffd@0:     N = size(dj,1);
wolffd@0: elseif size(dj,1)<N
wolffd@0:     dj(N,1,1) = 0;
wolffd@0: end
wolffd@0: if not(win == 0)
wolffd@0:     winf = str2func(win);
wolffd@0:     try
wolffd@0:         w = window(winf,N);
wolffd@0:     catch
wolffd@0:         if strcmpi(win,'hamming')
wolffd@0:             disp('Signal Processing Toolbox does not seem to be installed. Recompute the hamming window manually.');
wolffd@0:             w = 0.54 - 0.46 * cos(2*pi*(0:N-1)'/(N-1));
wolffd@0:         else
wolffd@0:             error(['ERROR in MIRSPECTRUM: Unknown windowing function ',win,' (maybe Signal Processing Toolbox is not installed).']);
wolffd@0:         end
wolffd@0:     end
wolffd@0:     kw = repmat(w,[1,size(dj,2),size(dj,3)]);
wolffd@0:     dj = dj(1:N,:,:).* kw;
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: function [y orig] = eachchunk(orig,option,missing,postchunk)
wolffd@0: option.zp = option.zp+missing;
wolffd@0: y = mirspectrum(orig,option);
wolffd@0: 
wolffd@0: 
wolffd@0: function y = combinechunk(old,new)
wolffd@0: do = get(old,'Data');
wolffd@0: do = do{1}{1};
wolffd@0: dn = get(new,'Data');
wolffd@0: dn = dn{1}{1};
wolffd@0: y = set(old,'ChunkData',do+dn);