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