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