Mercurial > hg > camir-aes2014

annotate toolboxes/MIRtoolbox1.3.2/MIRToolbox/@mirspectrum/mirspectrum.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 = mirspectrum(orig,varargin)
wolffd@0 2 % s = mirspectrum(x) computes the spectrum of the audio signal x, showing
wolffd@0 3 % the distribution of the energy along the frequencies.
wolffd@0 4 % (x can be the name of an audio file as well.)
wolffd@0 5 % Optional argument:
wolffd@0 6 % mirspectrum(...,'Frame',l,h) computes spectrogram, using frames of
wolffd@0 7 % length l seconds and a hop rate h.
wolffd@0 8 % Default values: l = .05 s, h = .5.
wolffd@0 9 % mirspectrum(...,'Min',mi) indicates the lowest frequency taken into
wolffd@0 10 % consideration, expressed in Hz.
wolffd@0 11 % Default value: 0 Hz.
wolffd@0 12 % mirspectrum(...,'Max',ma) indicates the highest frequency taken into
wolffd@0 13 % consideration, expressed in Hz.
wolffd@0 14 % Default value: the maximal possible frequency, corresponding to
wolffd@0 15 % the sampling rate of x divided by 2.
wolffd@0 16 % mirspectrum(...,'Window',w): windowing method:
wolffd@0 17 % either w = 0 (no windowing) or any windowing function proposed
wolffd@0 18 % in the Signal Processing Toolbox (help window).
wolffd@0 19 % default value: w = 'hamming'
wolffd@0 20 %
wolffd@0 21 % mirspectrum(...,'Cents'): decomposes the energy in cents.
wolffd@0 22 % mirspectrum(...,'Collapsed'): collapses the spectrum into one
wolffd@0 23 % octave divided into 1200 cents.
wolffd@0 24 % Redistribution of the frequencies into bands:
wolffd@0 25 % mirspectrum(...,'Mel'): Mel bands.
wolffd@0 26 % (Requires the Auditory Toolbox.)
wolffd@0 27 % mirspectrum(...,'Bark'): Bark bands.
wolffd@0 28 % (Code based on Pampalk's MA toolbox).
wolffd@0 29 % If the audio signal was frame decomposed, the output s is a
wolffd@0 30 % band-decomposed spectrogram. It is then possible to compute
wolffd@0 31 % the spectrum of the temporal signal in each band,
wolffd@0 32 % using the following syntax:
wolffd@0 33 % mirspectrum(s,'AlongBands')
wolffd@0 34 % This corresponds to fluctuation (cf. mirfluctuation).
wolffd@0 35 % mirspectrum(...,'Mask'): Models masking phenomena in each band.
wolffd@0 36 % (Code based on Pampalk's MA toolbox).
wolffd@0 37 % mirspectrum(...,'Normal'): normalizes with respect to energy.
wolffd@0 38 % mirspectrum(...,'NormalInput'): input signal is normalized from 0 to 1.
wolffd@0 39 % mirspectrum(...,'Power'): squares the energy.
wolffd@0 40 % mirspectrum(...,'dB'): represents the spectrum energy in decibel scale.
wolffd@0 41 % mirspectrum(...,'dB',th): keeps only the highest energy over a
wolffd@0 42 % range of th dB.
wolffd@0 43 % mirspectrum(...,'Terhardt'): modulates the energy following
wolffd@0 44 % (Terhardt, 1979) outer ear model. (Code based on Pampalk's MA
wolffd@0 45 % toolbox).
wolffd@0 46 % mirspectrum(...,'Resonance',r): multiplies the spectrum with a
wolffd@0 47 % resonance curve. Possible value for r:
wolffd@0 48 % 'ToiviainenSnyder': highlights best perceived meter
wolffd@0 49 % (Toiviainen & Snyder 2003)
wolffd@0 50 % (default value for spectrum of envelopes)
wolffd@0 51 % 'Fluctuation': fluctuation strength (Fastl 1982)
wolffd@0 52 % (default value for spectrum of band channels)
wolffd@0 53 % mirspectrum(...,'Prod',m): Enhances components that have harmonics
wolffd@0 54 % located at multiples of range(s) m of the signal's fundamental
wolffd@0 55 % frequency. Computed by compressing the signal by the list of
wolffd@0 56 % factors m, and by multiplying all the results with the original
wolffd@0 57 % signa (Alonso et al, 2003).
wolffd@0 58 % mirspectrum(...,'Sum',s): Similar option using summation instead of
wolffd@0 59 % product.
wolffd@0 60 %
wolffd@0 61 % mirspectrum(...,'MinRes',mr): Indicates a minimal accepted
wolffd@0 62 % frequency resolution, in Hz. The audio signal is zero-padded in
wolffd@0 63 % order to reach the desired resolution.
wolffd@0 64 % If the 'Mel' option is toggled on, 'MinRes' is set by default
wolffd@0 65 % to 66 Hz.
wolffd@0 66 % mirspectrum(...,'MinRes',mr,'OctaveRatio',tol): Indicates the
wolffd@0 67 % minimal accepted resolution in terms of number of divisions of
wolffd@0 68 % the octave. Low frequencies are ignored in order to reach the
wolffd@0 69 % desired resolution.
wolffd@0 70 % The corresponding required frequency resolution is equal to
wolffd@0 71 % the difference between the first frequency bins, multiplied
wolffd@0 72 % by the constraining multiplicative factor tol (set by
wolffd@0 73 % default to .75).
wolffd@0 74 % mirspectrum(...,'Res',r): Indicates the required precise frequency
wolffd@0 75 % resolution, in Hz. The audio signal is zero-padded in order to
wolffd@0 76 % reach the desired resolution.
wolffd@0 77 % mirspectrum(...,'Length',l): Specifies the length of the FFT,
wolffd@0 78 % overriding the FFT length initially planned.
wolffd@0 79 % mirspectrum(...,'ZeroPad',s): Zero-padding of s samples.
wolffd@0 80 % mirspectrum(...,'WarningRes',s): Indicates a required frequency
wolffd@0 81 % resolution, in Hz, for the input signal. If the resolution does
wolffd@0 82 % not reach that prerequisite, a warning is displayed.
wolffd@0 83 % mirspectrum(...,'ConstantQ',nb): Carries out a Constant Q Transform
wolffd@0 84 % instead of a FFT, with a number of bins per octave fixed to nb.
wolffd@0 85 % Default value for nb: 12 bins per octave.
wolffd@0 86 %
wolffd@0 87 % mirspectrum(...,'Smooth',o): smooths the envelope using a movering
wolffd@0 88 % average of order o.
wolffd@0 89 % Default value when the option is toggled on: o=10
wolffd@0 90 % mirspectrum(...,'Gauss',o): smooths the envelope using a gaussian
wolffd@0 91 % of standard deviation o samples.
wolffd@0 92 % Default value when the option is toggled on: o=10
wolffd@0 93 % mirspectrum(...,'Phase',0): do not compute the FFT phase.
wolffd@0 94
wolffd@0 95 win.key = 'Window';
wolffd@0 96 win.type = 'String';
wolffd@0 97 win.default = 'hamming';
wolffd@0 98 option.win = win;
wolffd@0 99
wolffd@0 100 min.key = 'Min';
wolffd@0 101 min.type = 'Integer';
wolffd@0 102 min.default = 0;
wolffd@0 103 option.min = min;
wolffd@0 104
wolffd@0 105 max.key = 'Max';
wolffd@0 106 max.type = 'Integer';
wolffd@0 107 max.default = Inf;
wolffd@0 108 option.max = max;
wolffd@0 109
wolffd@0 110 mr.key = 'MinRes';
wolffd@0 111 mr.type = 'Integer';
wolffd@0 112 mr.default = 0;
wolffd@0 113 option.mr = mr;
wolffd@0 114
wolffd@0 115 res.key = 'Res';
wolffd@0 116 res.type = 'Integer';
wolffd@0 117 res.default = NaN;
wolffd@0 118 option.res = res;
wolffd@0 119
wolffd@0 120 length.key = 'Length';
wolffd@0 121 length.type = 'Integer';
wolffd@0 122 length.default = NaN;
wolffd@0 123 option.length = length;
wolffd@0 124
wolffd@0 125 zp.key = 'ZeroPad';
wolffd@0 126 zp.type = 'Integer';
wolffd@0 127 zp.default = 0;
wolffd@0 128 zp.keydefault = Inf;
wolffd@0 129 option.zp = zp;
wolffd@0 130
wolffd@0 131 wr.key = 'WarningRes';
wolffd@0 132 wr.type = 'Integer';
wolffd@0 133 wr.default = 0;
wolffd@0 134 option.wr = wr;
wolffd@0 135
wolffd@0 136 octave.key = 'OctaveRatio';
wolffd@0 137 octave.type = 'Boolean';
wolffd@0 138 octave.default = 0;
wolffd@0 139 option.octave = octave;
wolffd@0 140
wolffd@0 141 constq.key = 'ConstantQ';
wolffd@0 142 constq.type = 'Integer';
wolffd@0 143 constq.default = 0;
wolffd@0 144 constq.keydefault = 12;
wolffd@0 145 option.constq = constq;
wolffd@0 146
wolffd@0 147 alongbands.key = 'AlongBands';
wolffd@0 148 alongbands.type = 'Boolean';
wolffd@0 149 alongbands.default = 0;
wolffd@0 150 option.alongbands = alongbands;
wolffd@0 151
wolffd@0 152 ni.key = 'NormalInput';
wolffd@0 153 ni.type = 'Boolean';
wolffd@0 154 ni.default = 0;
wolffd@0 155 option.ni = ni;
wolffd@0 156
wolffd@0 157 norm.key = 'Normal';
wolffd@0 158 norm.type = 'Integer';
wolffd@0 159 norm.default = 0;
wolffd@0 160 norm.keydefault = 1;
wolffd@0 161 norm.when = 'After';
wolffd@0 162 option.norm = norm;
wolffd@0 163
wolffd@0 164 band.type = 'String';
wolffd@0 165 band.choice = {'Freq','Mel','Bark','Cents'};
wolffd@0 166 band.default = 'Freq';
wolffd@0 167 band.when = 'Both';
wolffd@0 168 option.band = band;
wolffd@0 169
wolffd@0 170 nbbands.key = 'Bands';
wolffd@0 171 nbbands.type = 'Integer';
wolffd@0 172 nbbands.default = 0;
wolffd@0 173 nbbands.when = 'After';
wolffd@0 174 option.nbbands = nbbands;
wolffd@0 175
wolffd@0 176 mprod.key = 'Prod';
wolffd@0 177 mprod.type = 'Integers';
wolffd@0 178 mprod.default = [];
wolffd@0 179 mprod.keydefault = 2:6;
wolffd@0 180 mprod.when = 'After';
wolffd@0 181 option.mprod = mprod;
wolffd@0 182
wolffd@0 183 msum.key = 'Sum';
wolffd@0 184 msum.type = 'Integers';
wolffd@0 185 msum.default = [];
wolffd@0 186 msum.keydefault = 2:6;
wolffd@0 187 msum.when = 'After';
wolffd@0 188 option.msum = msum;
wolffd@0 189
wolffd@0 190 reso.key = 'Resonance';
wolffd@0 191 reso.type = 'String';
wolffd@0 192 reso.choice = {'ToiviainenSnyder','Fluctuation','Meter',0,'no','off'};
wolffd@0 193 reso.default = 0;
wolffd@0 194 reso.when = 'After';
wolffd@0 195 option.reso = reso;
wolffd@0 196
wolffd@0 197 log.key = 'log';
wolffd@0 198 log.type = 'Boolean';
wolffd@0 199 log.default = 0;
wolffd@0 200 log.when = 'After';
wolffd@0 201 option.log = log;
wolffd@0 202
wolffd@0 203 db.key = 'dB';
wolffd@0 204 db.type = 'Integer';
wolffd@0 205 db.default = 0;
wolffd@0 206 db.keydefault = Inf;
wolffd@0 207 db.when = 'After';
wolffd@0 208 option.db = db;
wolffd@0 209
wolffd@0 210 pow.key = 'Power';
wolffd@0 211 pow.type = 'Boolean';
wolffd@0 212 pow.default = 0;
wolffd@0 213 pow.when = 'After';
wolffd@0 214 option.pow = pow;
wolffd@0 215
wolffd@0 216 terhardt.key = 'Terhardt';
wolffd@0 217 terhardt.type = 'Boolean';
wolffd@0 218 terhardt.default = 0;
wolffd@0 219 terhardt.when = 'After';
wolffd@0 220 option.terhardt = terhardt;
wolffd@0 221
wolffd@0 222 mask.key = 'Mask';
wolffd@0 223 mask.type = 'Boolean';
wolffd@0 224 mask.default = 0;
wolffd@0 225 mask.when = 'After';
wolffd@0 226 option.mask = mask;
wolffd@0 227
wolffd@0 228 % e.key = 'Enhanced';
wolffd@0 229 % e.type = 'Integer';
wolffd@0 230 % e.default = [];
wolffd@0 231 % e.keydefault = 2:10;
wolffd@0 232 % e.when = 'After';
wolffd@0 233 %option.e = e;
wolffd@0 234
wolffd@0 235 collapsed.key = 'Collapsed';
wolffd@0 236 collapsed.type = 'Boolean';
wolffd@0 237 collapsed.default = 0;
wolffd@0 238 collapsed.when = 'Both';
wolffd@0 239 option.collapsed = collapsed;
wolffd@0 240
wolffd@0 241 aver.key = 'Smooth';
wolffd@0 242 aver.type = 'Integer';
wolffd@0 243 aver.default = 0;
wolffd@0 244 aver.keydefault = 10;
wolffd@0 245 aver.when = 'After';
wolffd@0 246 option.aver = aver;
wolffd@0 247
wolffd@0 248 gauss.key = 'Gauss';
wolffd@0 249 gauss.type = 'Integer';
wolffd@0 250 gauss.default = 0;
wolffd@0 251 gauss.keydefault = 10;
wolffd@0 252 gauss.when = 'After';
wolffd@0 253 option.gauss = gauss;
wolffd@0 254
wolffd@0 255 rapid.key = 'Rapid';
wolffd@0 256 rapid.type = 'Boolean';
wolffd@0 257 rapid.default = 0;
wolffd@0 258 option.rapid = rapid;
wolffd@0 259
wolffd@0 260 phase.key = 'Phase';
wolffd@0 261 phase.type = 'Boolean';
wolffd@0 262 phase.default = 1;
wolffd@0 263 option.phase = phase;
wolffd@0 264
wolffd@0 265 specif.option = option;
wolffd@0 266
wolffd@0 267 specif.defaultframelength = 0.05;
wolffd@0 268 specif.defaultframehop = 0.5;
wolffd@0 269
wolffd@0 270 specif.eachchunk = @eachchunk;
wolffd@0 271 specif.combinechunk = @combinechunk;
wolffd@0 272
wolffd@0 273 if isamir(orig,'mirscalar') || isamir(orig,'mirenvelope')
wolffd@0 274 specif.nochunk = 1;
wolffd@0 275 end
wolffd@0 276
wolffd@0 277 varargout = mirfunction(@mirspectrum,orig,varargin,nargout,specif,@init,@main);
wolffd@0 278
wolffd@0 279
wolffd@0 280 function [x type] = init(x,option)
wolffd@0 281 type = 'mirspectrum';
wolffd@0 282
wolffd@0 283
wolffd@0 284 function s = main(orig,option,postoption)
wolffd@0 285 if isstruct(option)
wolffd@0 286 if option.collapsed
wolffd@0 287 option.band = 'Cents';
wolffd@0 288 end
wolffd@0 289 if isnan(option.res) && strcmpi(option.band,'Cents') && option.min
wolffd@0 290 option.res = option.min *(2^(1/1200)-1)*.9;
wolffd@0 291 end
wolffd@0 292 end
wolffd@0 293 if not(isempty(postoption))
wolffd@0 294 if not(strcmpi(postoption.band,'Freq') && isempty(postoption.msum) ...
wolffd@0 295 || isempty(postoption.mprod)) || postoption.log || postoption.db ...
wolffd@0 296 || postoption.pow || postoption.mask || postoption.collapsed ...
wolffd@0 297 || postoption.aver || postoption.gauss
wolffd@0 298 option.phase = 0;
wolffd@0 299 end
wolffd@0 300 end
wolffd@0 301 if iscell(orig)
wolffd@0 302 orig = orig{1};
wolffd@0 303 end
wolffd@0 304 if isa(orig,'mirspectrum') && ...
wolffd@0 305 not(isfield(option,'alongbands') && option.alongbands)
wolffd@0 306 s = orig;
wolffd@0 307 if isfield(option,'min') && ...
wolffd@0 308 (option.min || iscell(option.max) || option.max < Inf)
wolffd@0 309 magn = get(s,'Magnitude');
wolffd@0 310 freq = get(s,'Frequency');
wolffd@0 311 for k = 1:length(magn)
wolffd@0 312 m = magn{k};
wolffd@0 313 f = freq{k};
wolffd@0 314 if iscell(option.max)
wolffd@0 315 mi = option.min{k};
wolffd@0 316 ma = option.max{k};
wolffd@0 317 else
wolffd@0 318 mi = option.min;
wolffd@0 319 ma = option.max;
wolffd@0 320 end
wolffd@0 321 if not(iscell(m))
wolffd@0 322 m = {m};
wolffd@0 323 f = {f};
wolffd@0 324 end
wolffd@0 325 for l = 1:length(m)
wolffd@0 326 range = find(and(f{l}(:,1) >= mi,f{l}(:,1) <= ma));
wolffd@0 327 m{l} = m{l}(range,:,:);
wolffd@0 328 f{l} = f{l}(range,:,:);
wolffd@0 329 end
wolffd@0 330 magn{k} = m;
wolffd@0 331 freq{k} = f;
wolffd@0 332 end
wolffd@0 333 s = set(s,'Magnitude',magn,'Frequency',freq);
wolffd@0 334 end
wolffd@0 335 if not(isempty(postoption)) && not(isequal(postoption,0))
wolffd@0 336 s = post(s,postoption);
wolffd@0 337 end
wolffd@0 338 elseif ischar(orig)
wolffd@0 339 s = mirspectrum(miraudio(orig),option,postoption);
wolffd@0 340 else
wolffd@0 341 if nargin == 0
wolffd@0 342 orig = [];
wolffd@0 343 end
wolffd@0 344 s.phase = [];
wolffd@0 345 s.log = 0;
wolffd@0 346 s.xscale = 'Freq';
wolffd@0 347 s.pow = 1;
wolffd@0 348 s = class(s,'mirspectrum',mirdata(orig));
wolffd@0 349 s = purgedata(s);
wolffd@0 350 s = set(s,'Title','Spectrum','Abs','frequency (Hz)','Ord','magnitude');
wolffd@0 351 %disp('Computing spectrum...')
wolffd@0 352 sig = get(orig,'Data');
wolffd@0 353 t = get(orig,'Pos');
wolffd@0 354 if isempty(t)
wolffd@0 355 t = get(orig,'FramePos');
wolffd@0 356 for k = 1:length(sig)
wolffd@0 357 for l = 1:length(sig{k})
wolffd@0 358 sig{k}{l} = sig{k}{l}';
wolffd@0 359 t{k}{l} = t{k}{l}(1,:,:)';
wolffd@0 360 end
wolffd@0 361 end
wolffd@0 362 end
wolffd@0 363 fs = get(orig,'Sampling');
wolffd@0 364 fp = get(orig,'FramePos');
wolffd@0 365 m = cell(1,length(sig));
wolffd@0 366 p = cell(1,length(sig));
wolffd@0 367 f = cell(1,length(sig));
wolffd@0 368 for i = 1:length(sig)
wolffd@0 369 d = sig{i};
wolffd@0 370 fpi = fp{i};
wolffd@0 371 if not(iscell(d))
wolffd@0 372 d = {d};
wolffd@0 373 end
wolffd@0 374 if option.alongbands
wolffd@0 375 fsi = 1 / (fpi{1}(1,2) - fpi{1}(1,1));
wolffd@0 376 else
wolffd@0 377 fsi = fs{i};
wolffd@0 378 end
wolffd@0 379 mi = cell(1,length(d));
wolffd@0 380 phi = cell(1,length(d));
wolffd@0 381 fi = cell(1,length(d));
wolffd@0 382 for J = 1:length(d)
wolffd@0 383 dj = d{J};
wolffd@0 384 if option.ni
wolffd@0 385 mxdj = repmat(max(dj),[size(dj,1),1,1]);
wolffd@0 386 mndj = repmat(min(dj),[size(dj,1),1,1]);
wolffd@0 387 dj = (dj-mndj)./(mxdj-mndj);
wolffd@0 388 end
wolffd@0 389 if option.alongbands
wolffd@0 390 if size(dj,1)>1
wolffd@0 391 error('ERROR IN MIRSPECTRUM: ''AlongBands'' is restricted to spectrum decomposed into bands, such as ''Mel'' and ''Bark''.')
wolffd@0 392 end
wolffd@0 393 dj = reshape(dj,[size(dj,2),1,size(dj,3)]);
wolffd@0 394 fp{i}{J} = fp{i}{J}([1;end]);
wolffd@0 395 end
wolffd@0 396
wolffd@0 397 if option.constq
wolffd@0 398 % Constant Q Transform
wolffd@0 399 r = 2^(1/option.constq);
wolffd@0 400 Q = 1 / (r - 1);
wolffd@0 401 f_max = min(fsi/2,option.max);
wolffd@0 402 f_min = option.min;
wolffd@0 403 if not(f_min)
wolffd@0 404 f_min = 16.3516;
wolffd@0 405 end
wolffd@0 406 B = floor(log(f_max/f_min) / log(r)); % number of bins
wolffd@0 407 N0 = round(Q*fsi/f_min); % maximum Nkcq
wolffd@0 408 j2piQn = -j*2*pi*Q*(0:N0-1)';
wolffd@0 409
wolffd@0 410 fj = f_min * r.^(0:B-1)';
wolffd@0 411 transf = NaN(B,size(dj,2),size(dj,3));
wolffd@0 412 for kcq = 1:B
wolffd@0 413 Nkcq = round(Q*fsi/fj(kcq));
wolffd@0 414 win = mirwindow(dj,option.win,Nkcq);
wolffd@0 415 exq = repmat(exp(j2piQn(1:Nkcq)/Nkcq),...
wolffd@0 416 [1,size(win,2),size(win,3)]);
wolffd@0 417 transf(kcq,:) = sum(win.* exq) / Nkcq;
wolffd@0 418 end
wolffd@0 419 else
wolffd@0 420 % FFT
wolffd@0 421 dj = mirwindow(dj,option.win);
wolffd@0 422 if option.zp
wolffd@0 423 if option.zp < Inf
wolffd@0 424 dj = [dj;zeros(option.zp,size(dj,2),size(dj,3))];
wolffd@0 425 else
wolffd@0 426 dj = [dj;zeros(size(dj))];
wolffd@0 427 end
wolffd@0 428 end
wolffd@0 429 if isstruct(postoption)
wolffd@0 430 if strcmpi(postoption.band,'Mel') && ...
wolffd@0 431 (not(option.mr) || option.mr > 66)
wolffd@0 432 option.mr = 66;
wolffd@0 433 end
wolffd@0 434 else
wolffd@0 435 %warning('WARNING in MIRSPECTRUM (for debugging purposes): By default, minimum resolution specified.')
wolffd@0 436 if not(option.mr)
wolffd@0 437 option.mr = 66;
wolffd@0 438 end
wolffd@0 439 end
wolffd@0 440 if option.octave
wolffd@0 441 N = size(dj,1);
wolffd@0 442 res = (2.^(1/option.mr)-1)*option.octave;
wolffd@0 443 % Minimal freq ratio between 2 first bins.
wolffd@0 444 % freq resolution should be > option.min * res
wolffd@0 445 Nrec = fsi/(option.min*res);
wolffd@0 446 % Recommended minimal sample length.
wolffd@0 447 if Nrec > N
wolffd@0 448 % If actual sample length is too small.
wolffd@0 449 option.min = fsi/N / res;
wolffd@0 450 warning('WARNING IN MIRSPECTRUM: The input signal is too short to obtain the desired octave resolution. Lowest frequencies will be ignored.');
wolffd@0 451 display(['(The recommended minimal input signal length would be ' num2str(Nrec/fsi) ' s.)']);
wolffd@0 452 display(['New low frequency range: ' num2str(option.min) ' Hz.']);
wolffd@0 453 end
wolffd@0 454 N = 2^nextpow2(N);
wolffd@0 455 elseif isnan(option.length)
wolffd@0 456 if isnan(option.res)
wolffd@0 457 N = size(dj,1);
wolffd@0 458 if option.mr && N < fsi/option.mr
wolffd@0 459 if option.wr && N < fsi/option.wr
wolffd@0 460 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 461 end
wolffd@0 462 N = max(N,fsi/option.mr);
wolffd@0 463 end
wolffd@0 464 N = 2^nextpow2(N);
wolffd@0 465 else
wolffd@0 466 N = ceil(fsi/option.res);
wolffd@0 467 end
wolffd@0 468 else
wolffd@0 469 N = option.length;
wolffd@0 470 end
wolffd@0 471
wolffd@0 472 % Here is the spectrum computation itself
wolffd@0 473 transf = fft(dj,N);
wolffd@0 474
wolffd@0 475 len = ceil(size(transf,1)/2);
wolffd@0 476 fj = fsi/2 * linspace(0,1,len)';
wolffd@0 477 if option.max
wolffd@0 478 maxf = find(fj>=option.max,1);
wolffd@0 479 if isempty(maxf)
wolffd@0 480 maxf = len;
wolffd@0 481 end
wolffd@0 482 else
wolffd@0 483 maxf = len;
wolffd@0 484 end
wolffd@0 485 if option.min
wolffd@0 486 minf = find(fj>=option.min,1);
wolffd@0 487 if isempty(minf)
wolffd@0 488 maxf = len;
wolffd@0 489 end
wolffd@0 490 else
wolffd@0 491 minf = 1;
wolffd@0 492 end
wolffd@0 493
wolffd@0 494 transf = transf(minf:maxf,:,:);
wolffd@0 495 fj = fj(minf:maxf);
wolffd@0 496 end
wolffd@0 497
wolffd@0 498 mi{J} = abs(transf);
wolffd@0 499 if option.phase
wolffd@0 500 phi{J} = angle(transf);
wolffd@0 501 end
wolffd@0 502 fi{J} = repmat(fj,[1,size(transf,2),size(transf,3)]);
wolffd@0 503 end
wolffd@0 504 if iscell(sig{i})
wolffd@0 505 m{i} = mi;
wolffd@0 506 p{i} = phi;
wolffd@0 507 f{i} = fi;
wolffd@0 508 else
wolffd@0 509 m{i} = mi{1};
wolffd@0 510 p{i} = phi{1};
wolffd@0 511 f{i} = fi{1};
wolffd@0 512 end
wolffd@0 513 end
wolffd@0 514 s = set(s,'Frequency',f,'Magnitude',m,'Phase',p,'FramePos',fp);
wolffd@0 515 if not(isempty(postoption)) && isstruct(postoption)
wolffd@0 516 s = post(s,postoption);
wolffd@0 517 end
wolffd@0 518 end
wolffd@0 519
wolffd@0 520
wolffd@0 521 function s = post(s,option)
wolffd@0 522 if option.collapsed
wolffd@0 523 option.band = 'Cents';
wolffd@0 524 end
wolffd@0 525 m = get(s,'Magnitude');
wolffd@0 526 f = get(s,'Frequency');
wolffd@0 527 for k = 1:length(m)
wolffd@0 528 if not(iscell(m{k}))
wolffd@0 529 m{k} = {m{k}};
wolffd@0 530 f{k} = {f{k}};
wolffd@0 531 end
wolffd@0 532 end
wolffd@0 533 if get(s,'Power') == 1 && ...
wolffd@0 534 (option.pow || any(option.mprod) || any(option.msum))
wolffd@0 535 for h = 1:length(m)
wolffd@0 536 for l = 1:length(m{k})
wolffd@0 537 m{h}{l} = m{h}{l}.^2;
wolffd@0 538 end
wolffd@0 539 end
wolffd@0 540 s = set(s,'Power',2,'Title',['Power ',get(s,'Title')]);
wolffd@0 541 end
wolffd@0 542 if any(option.mprod)
wolffd@0 543 for h = 1:length(m)
wolffd@0 544 for l = 1:length(m{k})
wolffd@0 545 z0 = m{h}{l};
wolffd@0 546 z1 = z0;
wolffd@0 547 for k = 1:length(option.mprod)
wolffd@0 548 mpr = option.mprod(k);
wolffd@0 549 if mpr
wolffd@0 550 zi = ones(size(z0));
wolffd@0 551 zi(1:floor(end/mpr),:,:) = z0(mpr:mpr:end,:,:);
wolffd@0 552 z1 = z1.*zi;
wolffd@0 553 end
wolffd@0 554 end
wolffd@0 555 m{h}{l} = z1;
wolffd@0 556 end
wolffd@0 557 end
wolffd@0 558 s = set(s,'Title','Spectral product');
wolffd@0 559 end
wolffd@0 560 if any(option.msum)
wolffd@0 561 for h = 1:length(m)
wolffd@0 562 for l = 1:length(m{k})
wolffd@0 563 z0 = m{h}{l};
wolffd@0 564 z1 = z0;
wolffd@0 565 for k = 1:length(option.msum)
wolffd@0 566 mpr = option.msum(k);
wolffd@0 567 if mpr
wolffd@0 568 zi = ones(size(z0));
wolffd@0 569 zi(1:floor(end/mpr),:,:) = z0(mpr:mpr:end,:,:);
wolffd@0 570 z1 = z1+zi;
wolffd@0 571 end
wolffd@0 572 end
wolffd@0 573 m{h}{l} = z1;
wolffd@0 574 end
wolffd@0 575 end
wolffd@0 576 s = set(s,'Title','Spectral sum');
wolffd@0 577 end
wolffd@0 578 %for i = option.e % Enhancing procedure %%%% TO CHECK, t IS NOT DEFINED!
wolffd@0 579 % for h = 1:length(m)
wolffd@0 580 % for l = 1:length(m{k})
wolffd@0 581 % c = m{h}{l};
wolffd@0 582 % % option.e is the list of scaling factors
wolffd@0 583 % % i is the scaling factor
wolffd@0 584 % if i
wolffd@0 585 % ic = zeros(size(c));
wolffd@0 586 % for g = 1:size(c,2)
wolffd@0 587 % for h2 = 1:size(c,3)
wolffd@0 588 % beg = find(t(:,g,h2)/i >= t(1,g,h2))
wolffd@0 589 % if not(isempty(beg))
wolffd@0 590 % ic(:,g,h2) = interp1(t(:,g,h2),c(:,g,h2),t(:,g,h2)/i); %,'cubic');
wolffd@0 591 % The scaled autocorrelation
wolffd@0 592 % ic(1:beg(1)-1,g,h2) = 0;
wolffd@0 593 % end
wolffd@0 594 % end
wolffd@0 595 % end
wolffd@0 596 % ic(find(isnan(ic))) = Inf; % All the NaN values are changed into 0 in the resulting curve
wolffd@0 597 % ic = max(ic,0);
wolffd@0 598 % c = c - ic; % The scaled autocorrelation is substracted to the initial one
wolffd@0 599 % c = max(c,0); % Half-wave rectification
wolffd@0 600 %figure
wolffd@0 601 %plot(c)
wolffd@0 602 % end
wolffd@0 603 % end
wolffd@0 604 % end
wolffd@0 605 %end
wolffd@0 606 if option.norm
wolffd@0 607 for k = 1:length(m)
wolffd@0 608 for l = 1:length(m{k})
wolffd@0 609 mkl = m{k}{l};
wolffd@0 610 nkl = zeros(1,size(mkl,2),size(mkl,3));
wolffd@0 611 for kk = 1:size(mkl,2)
wolffd@0 612 for ll = 1:size(mkl,3)
wolffd@0 613 nkl(1,kk,l) = norm(mkl(:,kk,ll));
wolffd@0 614 end
wolffd@0 615 end
wolffd@0 616 m{k}{l} = mkl./repmat(nkl,[size(m{k}{k},1),1,1]);
wolffd@0 617 end
wolffd@0 618 end
wolffd@0 619 end
wolffd@0 620 if option.terhardt && not(isempty(find(f{1}{1}))) % This excludes the case where spectrum already along bands
wolffd@0 621 % Code taken from Pampalk's MA Toolbox
wolffd@0 622 for h = 1:length(m)
wolffd@0 623 for l = 1:length(m{k})
wolffd@0 624 W_Adb = zeros(size(f{k}{l}));
wolffd@0 625 W_Adb(2:size(f{k}{l},1),:,:) = ...
wolffd@0 626 + 10.^((-3.64*(f{k}{l}(2:end,:,:)/1000).^-0.8 ...
wolffd@0 627 + 6.5 * exp(-0.6 * (f{k}{l}(2:end,:,:)/1000 - 3.3).^2) ...
wolffd@0 628 - 0.001*(f{k}{l}(2:end,:,:)/1000).^4)/20);
wolffd@0 629 W_Adb = W_Adb.^2;
wolffd@0 630 m{h}{l} = m{h}{l}.*W_Adb;
wolffd@0 631 end
wolffd@0 632 end
wolffd@0 633 end
wolffd@0 634 if option.reso
wolffd@0 635 if not(ischar(option.reso))
wolffd@0 636 if strcmp(get(orig,'XScale'),'Mel')
wolffd@0 637 option.reso = 'Fluctuation';
wolffd@0 638 else
wolffd@0 639 option.reso = 'ToiviainenSnyder';
wolffd@0 640 end
wolffd@0 641 end
wolffd@0 642 %if strcmpi(option.reso,'ToiviainenSnyder') || ...
wolffd@0 643 % strcmpi(option.reso,'Meter') || ...
wolffd@0 644 % strcmpi(option.reso,'Fluctuation')
wolffd@0 645 for k = 1:length(m)
wolffd@0 646 for l = 1:length(m{k})
wolffd@0 647 if strcmpi(option.reso,'ToiviainenSnyder') || strcmpi(option.reso,'Meter')
wolffd@0 648 w = max(0,...
wolffd@0 649 1 - 0.25*(log2(max(1./max(f{k}{l},1e-12),1e-12)/0.5)).^2);
wolffd@0 650 elseif strcmpi(option.reso,'Fluctuation')
wolffd@0 651 w1 = f{k}{l} / 4; % ascending part of the fluctuation curve;
wolffd@0 652 w2 = 1 - 0.3 * (f{k}{l} - 4)/6; % descending part;
wolffd@0 653 w = min(w1,w2);
wolffd@0 654 end
wolffd@0 655 if max(w) == 0
wolffd@0 656 warning('The resonance curve, not defined for this range of delays, will not be applied.')
wolffd@0 657 else
wolffd@0 658 m{k}{l} = m{k}{l}.* w;
wolffd@0 659 end
wolffd@0 660 end
wolffd@0 661 end
wolffd@0 662 end
wolffd@0 663 if strcmp(s.xscale,'Freq')
wolffd@0 664 if strcmpi(option.band,'Mel')
wolffd@0 665 % The following is largely based on the source code from Auditory Toolbox
wolffd@0 666 % (A part that I could not call directly from MIRtoolbox)
wolffd@0 667
wolffd@0 668 % (Malcolm Slaney, August 1993, (c) 1998 Interval Research Corporation)
wolffd@0 669
wolffd@0 670 try
wolffd@0 671 MakeERBFilters(1,1,1); % Just to be sure that the Auditory Toolbox is installed
wolffd@0 672 catch
wolffd@0 673 error('ERROR IN MIRFILTERBANK: Auditory Toolbox needs to be installed.');
wolffd@0 674 end
wolffd@0 675
wolffd@0 676 %disp('Computing Mel-frequency spectral representation...')
wolffd@0 677 lowestFrequency = 133.3333;
wolffd@0 678 if not(option.nbbands)
wolffd@0 679 option.nbbands = 40;
wolffd@0 680 end
wolffd@0 681 linearFilters = min(13,option.nbbands);
wolffd@0 682 linearSpacing = 66.66666666;
wolffd@0 683 logFilters = option.nbbands - linearFilters;
wolffd@0 684 logSpacing = 1.0711703;
wolffd@0 685 totalFilters = option.nbbands;
wolffd@0 686
wolffd@0 687 % Figure the band edges. Interesting frequencies are spaced
wolffd@0 688 % by linearSpacing for a while, then go logarithmic. First figure
wolffd@0 689 % all the interesting frequencies. Lower, center, and upper band
wolffd@0 690 % edges are all consequtive interesting frequencies.
wolffd@0 691 freqs = lowestFrequency + (0:linearFilters-1)*linearSpacing;
wolffd@0 692 freqs(linearFilters+1:totalFilters+2) = ...
wolffd@0 693 freqs(linearFilters) * logSpacing.^(1:logFilters+2);
wolffd@0 694 lower = freqs(1:totalFilters);
wolffd@0 695 center = freqs(2:totalFilters+1);
wolffd@0 696 upper = freqs(3:totalFilters+2);
wolffd@0 697
wolffd@0 698 % Figure out the height of the triangle so that each filter has
wolffd@0 699 % unit weight, assuming a triangular weighting function.
wolffd@0 700 triangleHeight = 2./(upper-lower);
wolffd@0 701
wolffd@0 702 e = cell(1,length(m));
wolffd@0 703 nch = cell(1,length(m));
wolffd@0 704 for h = 1:length(m)
wolffd@0 705 e{h} = cell(1,length(m{h}));
wolffd@0 706 for i = 1:length(m{h})
wolffd@0 707 mi = m{h}{i};
wolffd@0 708 fi = f{h}{i}(:,1,1);
wolffd@0 709 fftSize = size(fi,1);
wolffd@0 710
wolffd@0 711 % We now want to combine FFT bins and figure out
wolffd@0 712 % each frequencies contribution
wolffd@0 713 mfccFilterWeights = zeros(totalFilters,fftSize);
wolffd@0 714 for chan=1:totalFilters
wolffd@0 715 mfccFilterWeights(chan,:) = triangleHeight(chan).*...
wolffd@0 716 ((fi > lower(chan) & fi <= center(chan)).* ...
wolffd@0 717 (fi-lower(chan))/(center(chan)-lower(chan)) + ...
wolffd@0 718 (fi > center(chan) & fi < upper(chan)).* ...
wolffd@0 719 (upper(chan)-fi)/(upper(chan)-center(chan)));
wolffd@0 720 end
wolffd@0 721 if find(diff(not(sum(mfccFilterWeights,2)))==-1)
wolffd@0 722 % If one channel has no weight whereas the higher one
wolffd@0 723 % has a positive weight.
wolffd@0 724 warning('WARNING in MIRSPECTRUM: The frequency resolution of the spectrum is too low for the Mel transformation. Some Mel components are undefined.')
wolffd@0 725 display('Recommended frequency resolution: at least 66 Hz.')
wolffd@0 726 end
wolffd@0 727 e{h}{i} = zeros(1,size(mi,2),totalFilters);
wolffd@0 728 for J = 1:size(mi,2)
wolffd@0 729 if max(mi(:,J)) > 0
wolffd@0 730 fftData = zeros(fftSize,1); % Zero-padding ?
wolffd@0 731 fftData(1:size(mi,1)) = mi(:,J);
wolffd@0 732 p = mfccFilterWeights * fftData + 1e-16;
wolffd@0 733 e{h}{i}(1,J,:) = reshape(p,[1,1,length(p)]);
wolffd@0 734 end
wolffd@0 735 end
wolffd@0 736 f{h}{i} = zeros(1,size(mi,2),totalFilters);
wolffd@0 737 end
wolffd@0 738 nch{h} = (1:totalFilters)';
wolffd@0 739 end
wolffd@0 740 m = e;
wolffd@0 741 s = set(s,'XScale','Mel','Title','Mel-Spectrum','Abs','Mel bands','Channels',nch);
wolffd@0 742 elseif strcmpi(option.band,'Bark')
wolffd@0 743 sr = get(s,'Sampling');
wolffd@0 744 % Code taken from Pampalk's MA Toolbox.
wolffd@0 745 %% zwicker & fastl: psychoacoustics 1999, page 159
wolffd@0 746 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 747 e = cell(1,length(m));
wolffd@0 748 nch = cell(1,length(m));
wolffd@0 749 for h = 1:length(m)
wolffd@0 750 %% ignore critical bands outside of sampling rate range
wolffd@0 751 cb = min(min([find(bark_upper>sr{h}/2),length(bark_upper)]),length(bark_upper));
wolffd@0 752 e{h} = cell(1,length(m{h}));
wolffd@0 753 for i = 1:length(m{h})
wolffd@0 754 mi = sum(m{h}{i},3);
wolffd@0 755 e{h}{i} = zeros(1,size(mi,2),cb);
wolffd@0 756 k = 1;
wolffd@0 757 for J=1:cb, %% group into bark bands
wolffd@0 758 idx = find(f{h}{i}(k:end,1,1)<=bark_upper(J));
wolffd@0 759 idx = idx + k-1;
wolffd@0 760 e{h}{i}(1,:,J) = sum(mi(idx,:,:),1);
wolffd@0 761 k = max(idx)+1;
wolffd@0 762 end
wolffd@0 763 f{h}{i} = zeros(1,size(mi,2),cb);
wolffd@0 764 end
wolffd@0 765 nch{h} = (1:length(bark_upper))';
wolffd@0 766 end
wolffd@0 767 m = e;
wolffd@0 768 s = set(s,'XScale','Bark','Title','Bark-Spectrum','Abs','Bark bands','Channels',nch);
wolffd@0 769 elseif strcmpi(option.band,'Cents') || option.collapsed
wolffd@0 770 for h = 1:length(m)
wolffd@0 771 for i = 1:length(m{h})
wolffd@0 772 mi = m{h}{i};
wolffd@0 773 fi = f{h}{i};
wolffd@0 774 isgood = fi(:,1,1)*(2^(1/1200)-1) >= fi(2,1,1)-fi(1,1,1);
wolffd@0 775 good = find(isgood);
wolffd@0 776 if isempty(good)
wolffd@0 777 mirerror('mirspectrum',...
wolffd@0 778 'The frequency resolution of the spectrum is too low to be decomposed into cents.');
wolffd@0 779 display('Hint: if you specify a minimum value for the frequency range (''Min'' option)');
wolffd@0 780 display(' and if you do not specify any frequency resolution (''Res'' option), ');
wolffd@0 781 display(' the correct frequency resolution will be automatically chosen.');
wolffd@0 782 end
wolffd@0 783 if good>1
wolffd@0 784 warning(['MIRSPECTRUM: Frequency resolution in cent is achieved for frequencies over ',...
wolffd@0 785 num2str(floor(fi(good(1)))),' Hz. Lower frequencies are ignored.'])
wolffd@0 786 display('Hint: if you specify a minimum value for the frequency range (''Min'' option)');
wolffd@0 787 display(' and if you do not specify any frequency resolution (''Res'' option), ');
wolffd@0 788 display(' the correct frequency resolution will be automatically chosen.');
wolffd@0 789 end
wolffd@0 790 fi = fi(good,:,:);
wolffd@0 791 mi = mi(good,:,:);
wolffd@0 792 f2cents = 440*2.^(1/1200*(-1200*10:1200*10-1)');
wolffd@0 793 % The frequencies corresponding to the cent range
wolffd@0 794 cents = repmat((0:1199)',[20,size(fi,2),size(fi,3)]);
wolffd@0 795 % The cent range is for the moment centered to A440
wolffd@0 796 octaves = ones(1200,1)*(-10:10);
wolffd@0 797 octaves = repmat(octaves(:),[1,size(fi,2),size(fi,3)]);
wolffd@0 798 select = find(f2cents>fi(1) & f2cents<fi(end));
wolffd@0 799 % Cent range actually used in the current spectrum
wolffd@0 800 f2cents = repmat(f2cents(select),[1,size(fi,2),size(fi,3)]);
wolffd@0 801 cents = repmat(cents(select),[1,size(fi,2),size(fi,3)]);
wolffd@0 802 octaves = repmat(octaves(select),[1,size(fi,2),size(fi,3)]);
wolffd@0 803 m{h}{i} = interp1(fi(:,1,1),mi,f2cents(:,1,1));
wolffd@0 804 f{h}{i} = octaves*1200+cents + 6900;
wolffd@0 805 % Now the cent range is expressed in midicent
wolffd@0 806 end
wolffd@0 807 end
wolffd@0 808 s = set(s,'Abs','pitch (in midicents)','XScale','Cents');
wolffd@0 809 end
wolffd@0 810 end
wolffd@0 811 if option.mask
wolffd@0 812 if strcmp(s.xscale,'Freq')
wolffd@0 813 warning('WARNING IN MIRSPECTRUM: ''Mask'' option available only for Mel-spectrum and Bark-spectrum.');
wolffd@0 814 disp('''Mask'' option ignored');
wolffd@0 815 else
wolffd@0 816 nch = get(s,'Channels');
wolffd@0 817 for h = 1:length(m)
wolffd@0 818 % Code taken from Pampalk's MA Toolbox.
wolffd@0 819 %% spreading function: schroeder et al., 1979, JASA, optimizing
wolffd@0 820 %% digital speech coders by exploiting masking properties of the human ear
wolffd@0 821 cb = length(nch{h}); % Number of bands.
wolffd@0 822 for i = 1:cb,
wolffd@0 823 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 824 end
wolffd@0 825 for i = 1:length(m{h})
wolffd@0 826 for J = 1:size(m{h}{i},2)
wolffd@0 827 mj = m{h}{i}(1,J,:);
wolffd@0 828 mj = spread(1:length(mj),1:length(mj))*mj(:);
wolffd@0 829 m{h}{i}(1,J,:) = reshape(mj,1,1,length(mj));
wolffd@0 830 end
wolffd@0 831 end
wolffd@0 832 end
wolffd@0 833 end
wolffd@0 834 end
wolffd@0 835 if option.collapsed
wolffd@0 836 for h = 1:length(m)
wolffd@0 837 for i = 1:length(m{h})
wolffd@0 838 mi = m{h}{i};
wolffd@0 839 fi = f{h}{i};
wolffd@0 840 centclass = rem(fi(:,1,1),1200);
wolffd@0 841 %neg = find(centclass<0);
wolffd@0 842 %centclass(neg) = 1200 + centclass(neg);
wolffd@0 843 m{h}{i} = NaN(1200,size(mi,2),size(mi,3));
wolffd@0 844 for k = 0:1199
wolffd@0 845 m{h}{i}(k+1,:,:) = sum(mi(find(centclass==k),:,:),1);
wolffd@0 846 end
wolffd@0 847 f{h}{i} = repmat((0:1199)',[1,size(fi,2),size(fi,3)]);
wolffd@0 848 end
wolffd@0 849 end
wolffd@0 850 s = set(s,'Abs','Cents','XScale','Cents(Collapsed)');
wolffd@0 851 end
wolffd@0 852 if option.log || option.db
wolffd@0 853 if not(isa(s,'mirspectrum') && s.log)
wolffd@0 854 if option.db
wolffd@0 855 s = set(s,'log',10);
wolffd@0 856 else
wolffd@0 857 s = set(s,'log',1);
wolffd@0 858 end
wolffd@0 859 for k = 1:length(m)
wolffd@0 860 if not(iscell(m{k}))
wolffd@0 861 m{k} = {m{k}};
wolffd@0 862 f{k} = {f{k}};
wolffd@0 863 end
wolffd@0 864 for l = 1:length(m{k})
wolffd@0 865 m{k}{l} = log10(m{k}{l}+1e-16);
wolffd@0 866 if option.db
wolffd@0 867 m{k}{l} = 10*m{k}{l};
wolffd@0 868 end
wolffd@0 869 end
wolffd@0 870 end
wolffd@0 871 elseif isa(s,'mirspectrum') && s.log<10
wolffd@0 872 for k = 1:length(m)
wolffd@0 873 for l = 1:length(m{k})
wolffd@0 874 m{k}{l} = 10*m{k}{l};
wolffd@0 875 end
wolffd@0 876 end
wolffd@0 877 end
wolffd@0 878 if option.db>0 && option.db < Inf
wolffd@0 879 for k = 1:length(m)
wolffd@0 880 for l = 1:length(m{k})
wolffd@0 881 m{k}{l} = m{k}{l}-repmat(max(m{k}{l}),[size(m{k}{l},1) 1 1]);
wolffd@0 882 m{k}{l} = option.db + max(-option.db,m{k}{l});
wolffd@0 883 end
wolffd@0 884 end
wolffd@0 885 end
wolffd@0 886 end
wolffd@0 887 if option.aver
wolffd@0 888 for k = 1:length(m)
wolffd@0 889 for i = 1:length(m{k})
wolffd@0 890 m{k}{i} = filter(ones(1,option.aver),1,m{k}{i});
wolffd@0 891 end
wolffd@0 892 end
wolffd@0 893 end
wolffd@0 894 if option.gauss
wolffd@0 895 for k = 1:length(m)
wolffd@0 896 for i = 1:length(m{k})
wolffd@0 897 sigma = option.gauss;
wolffd@0 898 gauss = 1/sigma/2/pi...
wolffd@0 899 *exp(- (-4*sigma:4*sigma).^2 /2/sigma^2);
wolffd@0 900 y = filter(gauss,1,[m{k}{i};zeros(4*sigma,1)]);
wolffd@0 901 y = y(4*sigma:end,:,:);
wolffd@0 902 m{k}{i} = y(1:size(m{k}{i},1),:,:);
wolffd@0 903 end
wolffd@0 904 end
wolffd@0 905 end
wolffd@0 906 s = set(s,'Magnitude',m,'Frequency',f);
wolffd@0 907
wolffd@0 908
wolffd@0 909 function dj = mirwindow(dj,win,N)
wolffd@0 910 if nargin<3
wolffd@0 911 N = size(dj,1);
wolffd@0 912 elseif size(dj,1)<N
wolffd@0 913 dj(N,1,1) = 0;
wolffd@0 914 end
wolffd@0 915 if not(win == 0)
wolffd@0 916 winf = str2func(win);
wolffd@0 917 try
wolffd@0 918 w = window(winf,N);
wolffd@0 919 catch
wolffd@0 920 if strcmpi(win,'hamming')
wolffd@0 921 disp('Signal Processing Toolbox does not seem to be installed. Recompute the hamming window manually.');
wolffd@0 922 w = 0.54 - 0.46 * cos(2*pi*(0:N-1)'/(N-1));
wolffd@0 923 else
wolffd@0 924 error(['ERROR in MIRSPECTRUM: Unknown windowing function ',win,' (maybe Signal Processing Toolbox is not installed).']);
wolffd@0 925 end
wolffd@0 926 end
wolffd@0 927 kw = repmat(w,[1,size(dj,2),size(dj,3)]);
wolffd@0 928 dj = dj(1:N,:,:).* kw;
wolffd@0 929 end
wolffd@0 930
wolffd@0 931
wolffd@0 932 function [y orig] = eachchunk(orig,option,missing,postchunk)
wolffd@0 933 option.zp = option.zp+missing;
wolffd@0 934 y = mirspectrum(orig,option);
wolffd@0 935
wolffd@0 936
wolffd@0 937 function y = combinechunk(old,new)
wolffd@0 938 do = get(old,'Data');
wolffd@0 939 do = do{1}{1};
wolffd@0 940 dn = get(new,'Data');
wolffd@0 941 dn = dn{1}{1};
wolffd@0 942 y = set(old,'ChunkData',do+dn);