v_spgrambw

PURPOSE ^

V_SPGRAMBW Draw spectrogram [T,F,B]=(s,fs,mode,bw,fmax,db,tinc,ann)

SYNOPSIS ^

function [t,f,b]=v_spgrambw(s,fs,varargin)

DESCRIPTION ^

V_SPGRAMBW Draw spectrogram [T,F,B]=(s,fs,mode,bw,fmax,db,tinc,ann)

  Usage: (1) v_spgrambw(s,fs,'pJcw')                       % Plot spectrogram with my favourite set of options

         (2) [y,fs,wrd,phn]=v_readsph(filename);           % read a TIMIT file
             v_spgrambw(s,fs,'pJcwat',[],[],[],[],wrd);    % plot spectrogram with transcription (replace
                                                             wrd by phn for phonetic trascription)

         (3) v_spgrambw(s,fs,'PJcwm',50,[100 2000])        % Plot narrow-band spectrogram on mel scale
                                                             from 100 to 2000 mel in power/mel units

         (4) [t,f,b]=v_spgrambw(s,fs,'p');                 % calculate the spectrogram without plotting
             imagesc(t,f,10*log10(b'));                    % plot it manually
             axis('xy');

         (5) ninc=0.0045*fs;                               % Frame increment for BW=200 Hz (in samples)
             nwin=2*ninc;                                  % Frame length (in samples)
             win=hamming(nwin);                            % Analysis window
             k=0.5*fs*sum(win.^2);                         % Scale factor to convert to power/Hz
             sf=abs(v_rfft(v_enframe(s,win,ninc),nwin,2)).^2/k;         % Calculate spectrum array                
             v_spgrambw(sf,[fs/ninc 0.5*(nwin+1)/fs fs/nwin],'Jc',bw);  % Plot spectrum array

         For examples of the many options available see:
         http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/tutorial/spgrambw/spgram_tut.pdf

  Inputs:  S         speech signal, or single-sided power spectrum array, S(NT,NF), in power per Hz
           FS        sample fequency (Hz) or [FS T1] where T1 is the time of the first sample
                     or, if s is a matrix, [FS T1 FINC F1] where FS is the frame rate, T1 is
                     the time of the first sample, FINC is the frequency increment and F1 the
                     frequency of the first column.
           MODE      optional character string specifying options (see list below)
           BW        bandwidth resolution in Hz (DFT window length = 1.81/BW)[default: 200]
           FMAX      frequency range [Fmin Fstep Fmax]. If Fstep is omitted
                     it is taken to be (Fmax-Fmin)/257, if Fmin is also omitted it is taken
                     to be 0 (or 20Hz for mode l), if all three are omitted Fmax is taken to be FS/2.
                     If modes m, b, e or l are specified then the units are in mel, bark or erb or
                     log10(Hz); this can be over-ridden by the 'h' option.
           DB        either dB-range relative to peak or [dB-min dB-max] for plotting (and B output if 'D' given [default: 40]
           TINC      output frame increment in seconds [0 or missing uses default=0.45/BW]
                     or [TFIRST TLAST] or [TFIRST TINC TLAST] where TFIRST/TLAST are the times
                     of first/last frames
           ANN       annotation cell array: each row contains either
                     {time 'text-string' 'font'} or {[t_start t_end] 'text-string' 'font'} where
                     the time value is in seconds with s(n) at time offset+n/fs. The font column can
                     omitted in which case the system font will be used. MATLAB cannot cope with
                     unicode so I recommend the SILDoulosIPA (serifed) or SILSophiaIPA (sans) fonts
                     for phonetic symbols; these are now a little hard to find.

 Outputs:  T(NT)        time axis values (in seconds). Input sample s(n) is at time offset+n/fs.
           F(NF)        frequency axis values in Hz or, unless mode=H, other selected frequency units
                        according to mode: m=mel, l=log10(Hz), b=bark,e=erb-rate
           B(NT,NF)     spectrogram values in power per x where x depends on the 'pPmbel' options
                        clipped to DB range if 'D' option and in dB if 'd' option.

 MODE:  'p' = output power per decade rather than power per Hz [preemphasis]
        'P' = output power per mel/bark/erb according to y axis scaling
        'd' = output B array is in dB rather than power
        'D' = clip the output B array to the limits specified by the "db" input

        'n' = use nearest frequency bin instead of interpolating

        'm' = mel scale
        'b' = bark scale
        'e' = erb scale
        'l' = log10 Hz frequency scale
        'f' = label frequency axis in Hz rather than mel/bark/...

        'h' = units of the FMAX input are in Hz instead of mel/bark
              [in this case, the Fstep parameter is used only to determine
               the number of filters]
        'H' = express the F output in Hz instead of mel/bark/...

        'g' = draw a graph even if output arguments are present
        'j' = jet colourmap
        'J' = "thermal" colourmap that is linear in grayscale. Based on Oliver Woodford's
                 real2rgb at http://www.mathworks.com/matlabcentral/fileexchange/23342
        'i' = inverted colourmap (white background)
        'c' = include a colourbar as an intensity scale
        'w' = draw the speech waveform above the spectrogram
        'a' = centre-align annotations rather than left-aligning them
        't' = add time markers with annotations

 The BW input gives the 6dB bandwidth of the Hamming window used in the analysis.
 Equal amplitude frequency components are guaranteed to give separate peaks if they
 are this far apart. This value also determines the time resolution: the window length is
 1.81/BW and the low-pass filter applied to amplitude modulations has a 6-dB bandwidth of
 BW/2 Hz.

 The units are power per Hz unless the 'P' option is given in which case power
 per displayed unit is used or power per decade for the 'l' option.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [t,f,b]=v_spgrambw(s,fs,varargin)
0002 %V_SPGRAMBW Draw spectrogram [T,F,B]=(s,fs,mode,bw,fmax,db,tinc,ann)
0003 %
0004 %  Usage: (1) v_spgrambw(s,fs,'pJcw')                       % Plot spectrogram with my favourite set of options
0005 %
0006 %         (2) [y,fs,wrd,phn]=v_readsph(filename);           % read a TIMIT file
0007 %             v_spgrambw(s,fs,'pJcwat',[],[],[],[],wrd);    % plot spectrogram with transcription (replace
0008 %                                                             wrd by phn for phonetic trascription)
0009 %
0010 %         (3) v_spgrambw(s,fs,'PJcwm',50,[100 2000])        % Plot narrow-band spectrogram on mel scale
0011 %                                                             from 100 to 2000 mel in power/mel units
0012 %
0013 %         (4) [t,f,b]=v_spgrambw(s,fs,'p');                 % calculate the spectrogram without plotting
0014 %             imagesc(t,f,10*log10(b'));                    % plot it manually
0015 %             axis('xy');
0016 %
0017 %         (5) ninc=0.0045*fs;                               % Frame increment for BW=200 Hz (in samples)
0018 %             nwin=2*ninc;                                  % Frame length (in samples)
0019 %             win=hamming(nwin);                            % Analysis window
0020 %             k=0.5*fs*sum(win.^2);                         % Scale factor to convert to power/Hz
0021 %             sf=abs(v_rfft(v_enframe(s,win,ninc),nwin,2)).^2/k;         % Calculate spectrum array
0022 %             v_spgrambw(sf,[fs/ninc 0.5*(nwin+1)/fs fs/nwin],'Jc',bw);  % Plot spectrum array
0023 %
0024 %         For examples of the many options available see:
0025 %         http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/tutorial/spgrambw/spgram_tut.pdf
0026 %
0027 %  Inputs:  S         speech signal, or single-sided power spectrum array, S(NT,NF), in power per Hz
0028 %           FS        sample fequency (Hz) or [FS T1] where T1 is the time of the first sample
0029 %                     or, if s is a matrix, [FS T1 FINC F1] where FS is the frame rate, T1 is
0030 %                     the time of the first sample, FINC is the frequency increment and F1 the
0031 %                     frequency of the first column.
0032 %           MODE      optional character string specifying options (see list below)
0033 %           BW        bandwidth resolution in Hz (DFT window length = 1.81/BW)[default: 200]
0034 %           FMAX      frequency range [Fmin Fstep Fmax]. If Fstep is omitted
0035 %                     it is taken to be (Fmax-Fmin)/257, if Fmin is also omitted it is taken
0036 %                     to be 0 (or 20Hz for mode l), if all three are omitted Fmax is taken to be FS/2.
0037 %                     If modes m, b, e or l are specified then the units are in mel, bark or erb or
0038 %                     log10(Hz); this can be over-ridden by the 'h' option.
0039 %           DB        either dB-range relative to peak or [dB-min dB-max] for plotting (and B output if 'D' given [default: 40]
0040 %           TINC      output frame increment in seconds [0 or missing uses default=0.45/BW]
0041 %                     or [TFIRST TLAST] or [TFIRST TINC TLAST] where TFIRST/TLAST are the times
0042 %                     of first/last frames
0043 %           ANN       annotation cell array: each row contains either
0044 %                     {time 'text-string' 'font'} or {[t_start t_end] 'text-string' 'font'} where
0045 %                     the time value is in seconds with s(n) at time offset+n/fs. The font column can
0046 %                     omitted in which case the system font will be used. MATLAB cannot cope with
0047 %                     unicode so I recommend the SILDoulosIPA (serifed) or SILSophiaIPA (sans) fonts
0048 %                     for phonetic symbols; these are now a little hard to find.
0049 %
0050 % Outputs:  T(NT)        time axis values (in seconds). Input sample s(n) is at time offset+n/fs.
0051 %           F(NF)        frequency axis values in Hz or, unless mode=H, other selected frequency units
0052 %                        according to mode: m=mel, l=log10(Hz), b=bark,e=erb-rate
0053 %           B(NT,NF)     spectrogram values in power per x where x depends on the 'pPmbel' options
0054 %                        clipped to DB range if 'D' option and in dB if 'd' option.
0055 %
0056 % MODE:  'p' = output power per decade rather than power per Hz [preemphasis]
0057 %        'P' = output power per mel/bark/erb according to y axis scaling
0058 %        'd' = output B array is in dB rather than power
0059 %        'D' = clip the output B array to the limits specified by the "db" input
0060 %
0061 %        'n' = use nearest frequency bin instead of interpolating
0062 %
0063 %        'm' = mel scale
0064 %        'b' = bark scale
0065 %        'e' = erb scale
0066 %        'l' = log10 Hz frequency scale
0067 %        'f' = label frequency axis in Hz rather than mel/bark/...
0068 %
0069 %        'h' = units of the FMAX input are in Hz instead of mel/bark
0070 %              [in this case, the Fstep parameter is used only to determine
0071 %               the number of filters]
0072 %        'H' = express the F output in Hz instead of mel/bark/...
0073 %
0074 %        'g' = draw a graph even if output arguments are present
0075 %        'j' = jet colourmap
0076 %        'J' = "thermal" colourmap that is linear in grayscale. Based on Oliver Woodford's
0077 %                 real2rgb at http://www.mathworks.com/matlabcentral/fileexchange/23342
0078 %        'i' = inverted colourmap (white background)
0079 %        'c' = include a colourbar as an intensity scale
0080 %        'w' = draw the speech waveform above the spectrogram
0081 %        'a' = centre-align annotations rather than left-aligning them
0082 %        't' = add time markers with annotations
0083 %
0084 % The BW input gives the 6dB bandwidth of the Hamming window used in the analysis.
0085 % Equal amplitude frequency components are guaranteed to give separate peaks if they
0086 % are this far apart. This value also determines the time resolution: the window length is
0087 % 1.81/BW and the low-pass filter applied to amplitude modulations has a 6-dB bandwidth of
0088 % BW/2 Hz.
0089 %
0090 % The units are power per Hz unless the 'P' option is given in which case power
0091 % per displayed unit is used or power per decade for the 'l' option.
0092 
0093 %%%% BUGS %%%%%%
0094 % * allow ANN rows to be a mixture of intervals and instants
0095 % * allow multiple ANN rows
0096 % * Do not use triangular interpolation if the output frequencies are the same as an FFT
0097 % * Place as many subticks as will fit beyond the last tick with the 'f' option
0098 % * Use a special subtick pattern between ticks that are powers of 10 using the 'f' option
0099 % * Future options:
0100 %       ['q' = constant q transform]
0101 %       ['k' = add a piano keyboard to the frequency scale]
0102 %       ['z' = use a bipolar colourmap for a matrix input with negative values]
0103 
0104 %      Copyright (C) Mike Brookes 1997-2011
0105 %      Version: $Id: v_spgrambw.m 10865 2018-09-21 17:22:45Z dmb $
0106 %
0107 %   VOICEBOX is a MATLAB toolbox for speech processing.
0108 %   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
0109 %
0110 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0111 %   This program is free software; you can redistribute it and/or modify
0112 %   it under the terms of the GNU General Public License as published by
0113 %   the Free Software Foundation; either version 2 of the License, or
0114 %   (at your option) any later version.
0115 %
0116 %   This program is distributed in the hope that it will be useful,
0117 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0118 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0119 %   GNU General Public License for more details.
0120 %
0121 %   You can obtain a copy of the GNU General Public License from
0122 %   http://www.gnu.org/copyleft/gpl.html or by writing to
0123 %   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
0124 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0125 persistent tcmap
0126 if isempty(tcmap)
0127     % modified thermal with better grayscale linearity
0128     tcmap=[ 0 0 0; 7 0 17; 14 0 33; 21 0 50; 29 0 67; 36 0 84; 43 0 100; 50 0 117;
0129         57 0 134; 64 0 150; 72 0 167; 80 3 164; 89 7 156; 97 11 149; 106 15 142; 114 19 134;
0130         123 23 127; 131 27 119; 140 31 112; 149 35 105; 157 39 97; 166 43 90; 174 47 82;
0131         183 51 75; 192 55 68; 200 59 60; 209 63 53; 217 67 45; 226 71 38; 234 75 31;
0132         243 79 23; 252 83 16; 255 88 12; 255 95 12; 255 102 11; 255 109 11; 255 116 10;
0133         255 123 10; 255 130 9; 255 137 9; 255 144 8; 255 151 8; 255 158 7; 255 165 7;
0134         255 172 6; 255 179 6; 255 186 5; 255 193 4; 255 200 4; 255 207 3; 255 214 3; 255 221 2;
0135         255 228 2; 255 235 1; 255 242 1; 255 249 0; 255 252 22; 255 252 55; 255 253 88;
0136         255 253 122; 255 254 155; 255 254 188; 255 255 222; 255 255 255]/255;
0137 end
0138 if nargin<2
0139     error('Usage: SPGRAMBW(s,fs,mode,bw,fmax,db,tinc)');
0140 end
0141 %SPGRAMBW Draw grey-scale spectrogram [T,F,B]=(s,fs,mode,bw,fmax,db,tinc)
0142 %
0143 % first decode the input arguments
0144 %
0145 if size(s,1)==1
0146     s=s(:);   % force to be a column vector (unless it is a matrix)
0147 end
0148 [ns1,ns2]=size(s);
0149 ap=zeros(1,6);
0150 j=2;
0151 if numel(fs)<2
0152     fs(2)=1/fs(1);  % first sample or frame is at time 1/fs
0153 end
0154 for i=1:length(varargin)
0155     if ischar(varargin{i})
0156         ap(1)=i;
0157     else
0158         ap(j)=i;
0159         j=j+1;
0160     end
0161 end
0162 if ap(1) && ~isempty(varargin{ap(1)})
0163     mode=varargin{ap(1)};
0164 else
0165     mode='';  % default mode
0166 end
0167 if ap(2) && ~isempty(varargin{ap(2)})
0168     bw=varargin{ap(2)};
0169 else
0170     bw=200;
0171 end
0172 if ap(3) && ~isempty(varargin{ap(3)})
0173     fmax=varargin{ap(3)};
0174 else
0175     fmax=[];
0176 end
0177 if ap(4) && ~isempty(varargin{ap(4)})
0178     db=varargin{ap(4)};
0179 else
0180     db=40;
0181 end
0182 if ap(5) && ~isempty(varargin{ap(5)})
0183     tinc=varargin{ap(5)};
0184 else
0185     tinc=0;
0186 end
0187 switch numel(tinc)
0188     case 1
0189         tinc=[tinc -Inf Inf];
0190     case 2
0191         tinc=[0 tinc];
0192     otherwise
0193         tinc=tinc([2 1 3]);
0194 end
0195 if tinc(1)<=0
0196     tinc(1)=1.81/(4*bw); % default frame increment
0197 end
0198 if ap(6)
0199     ann=varargin{ap(6)};
0200 else
0201     ann=[];
0202 end
0203 
0204 % now sort out the mode flags
0205 
0206 mdsw='  ';           % [yscale preemph]
0207 for i=1:length(mode)
0208     switch mode(i)
0209         case {'l','m','b','e'}
0210             mdsw(1)=mode(i);
0211         case {'p','P'}
0212             mdsw(2)=mode(i);
0213     end
0214 end
0215 if mdsw(2)=='P'
0216     mdsw(2)=mdsw(1);        % preemphasis is scaling dependent
0217 end
0218 %
0219 % sort out the frequency axis
0220 %
0221 flmin=30;                   % min frequency for 'l' option
0222 nfrq=257;                   % default number of frequency bins
0223 if ns2==1
0224     fnyq=fs(1)/2;           % default upper frequency limit is fs/2
0225 else                        % input is a power spectrum
0226     if numel(fs)<3
0227         fs(3)=fs(1)*0.25;   % default increment is 0.25 times frame increment
0228     end
0229     if numel(fs)<4
0230         fs(4)=0;            % first freq bin is DC by default
0231     end
0232     fnyq=fs(4)+(ns2-1)*fs(3);  % default upper frequency limit is highest supplied frequency
0233 end
0234 
0235 if ~numel(fmax)             % no explicit frequency range
0236     switch mdsw(1)
0237         case 'l'
0238             fx=linspace(log10(flmin),log10(fnyq),nfrq);   % 20  Hz to Nyquist
0239         case 'm'
0240             fx=linspace(0,v_frq2mel(fnyq),nfrq);   % DC to Nyquist
0241         case 'b'
0242             fx=linspace(0,v_frq2bark(fnyq),nfrq);   % DC to Nyquist
0243         case 'e'
0244             fx=linspace(0,v_frq2erb(fnyq),nfrq);   % DC to Nyquist
0245         otherwise   % linear Hz scale
0246             fx=(0:nfrq-1)*fnyq/(nfrq-1);
0247     end
0248 else
0249     if any(mode=='h')
0250         switch mdsw(1)
0251             case 'l'
0252                 fmaxu=log10(fmax);   % 20  Hz to Nyquist
0253             case 'm'
0254                 fmaxu=v_frq2mel(fmax);   % DC to Nyquist
0255             case 'b'
0256                 fmaxu=v_frq2bark(fmax);   % DC to Nyquist
0257             case 'e'
0258                 fmaxu=v_frq2erb(fmax);   % DC to Nyquist
0259             otherwise
0260                 fmaxu=fmax;  % linear Hz scale
0261         end
0262     else
0263         fmaxu=fmax;                 % already in the correct units
0264     end
0265     if numel(fmax)<2   % only max value specified
0266         if mdsw(1)=='l'
0267             fx=linspace(log10(flmin),fmaxu,nfrq);   % 20  Hz to fmax
0268         else
0269             fx=linspace(0,fmaxu,nfrq);   % DC to fmax
0270         end
0271     elseif numel(fmax)<3 % min and max values specified
0272         fx=linspace(fmaxu(1),fmaxu(2),nfrq);   % fmin to fmax
0273     else
0274         fmaxu(2)=fmax(2)*(fmaxu(3)-fmaxu(1))/(fmax(3)-fmax(1)); % scale the step size appropriately
0275         fx=fmaxu(1):fmaxu(2):fmaxu(3);   % fmin to fmax in steps of finc
0276         nfrq=length(fx);
0277     end
0278 end
0279 switch mdsw(1)          % convert the frequency range to Hz
0280     case 'l'
0281         f=10.^fx;
0282         frlab='log_{10}Hz';
0283         frlabf='log';
0284         frq2y=@log10;
0285         y2frq=@(x) 10.^x;
0286     case 'm'
0287         f=v_mel2frq(fx);
0288         frlab='Mel';
0289         frlabf='Mel';
0290         frq2y=@v_frq2mel;
0291         y2frq=@v_mel2frq;
0292     case 'b'
0293         f=v_bark2frq(fx);
0294         frlab='Bark';
0295         frlabf='Bark';
0296         frq2y=@v_frq2bark;
0297         y2frq=@v_bark2frq;
0298     case 'e'
0299         f=v_erb2frq(fx);
0300         frlab='Erb-rate';
0301         frlabf='Erb';
0302         frq2y=@v_frq2erb;
0303         y2frq=@v_erb2frq;
0304     otherwise
0305         f=fx;
0306         frlab='Hz';
0307         frq2y=@(x) x;
0308         y2frq=@(x) x;
0309 end
0310 if ~any(mode=='H')
0311     f=fx;               % give output frequencies in native units instead of Hz unless 'H' is specified
0312 end
0313 %
0314 % now calculate the spectrogram
0315 %
0316 if ns2==1   % input is a speech signal vector
0317     winlen = fix(1.81*fs(1)/bw);   % window length
0318     win=0.54+0.46*cos((1-winlen:2:winlen)*pi/winlen);  % Hamming window
0319     ninc=max(round(tinc(1)*fs(1)),1);                 % window increment in samples
0320     %  we need to take account of minimum freq increment + make it exact if possible
0321     fftlen=pow2(nextpow2(4*winlen));        % enough oversampling to get good interpolation
0322     win=win/sqrt(sum(win.^2));              % ensure window squared sums to unity
0323     ix1=max(round((tinc(2)-fs(2))*fs(1)-(winlen-3)/2),1); % first sample required
0324     ix2=min(ceil((tinc(3)-fs(2))*fs(1)+(winlen+1)/2),ns1); % last sample required
0325     [sf,t]=v_enframe(s(ix1:ix2),win,ninc);
0326     t=fs(2)+(t+ix1-2)/fs(1);                         % time axis
0327     b=v_rfft(sf,fftlen,2);
0328     b=b.*conj(b)*2/fs(1);          % Power per Hz
0329     b(:,1)=b(:,1)*0.5;   % correct for no negative zero frequency to double the power
0330     b(:,end)=b(:,end)*0.5;   % correct for no negative nyquist frequency to double the power
0331     fb=(0:fftlen/2)*fs(1)/fftlen; % fft bin frequencies
0332     fftfs=fs(1);
0333 else
0334     b=s;
0335     t=fs(2)+(0:ns1-1)/fs(1);  % frame times
0336     fb=fs(4)+(0:ns2-1)*fs(3);
0337     fftlen=[ns2 fs(3) fs(4)]; % for v_filtbankm: ns2=# input freq bins, freq increment (Hz), first bin freq (Hz)
0338     fftfs=0;
0339     %     fftlen=2*(ns2-1);  % assume an even length fft
0340     %     fftfs=fftlen*fs(3);
0341 end
0342 nfr=numel(t);                   % number of frames
0343 dblab='Power/Hz';
0344 switch mdsw(2)
0345     case {'p','l'}
0346         b=b.*repmat(fb*log(10),nfr,1);       % convert to power per decade
0347         dblab='Power/Decade';
0348     case 'm'
0349         b=b.*repmat((700+fb)*log(1+1000/700)/1000,nfr,1);       % convert to power per mel
0350         dblab='Power/Mel';
0351     case 'b'
0352         b=b.*repmat((1960+fb).^2/52547.6,nfr,1);       % convert to power per bark
0353         dblab='Power/Bark';
0354     case 'e'
0355         b=b.*repmat(6.23*fb.^2 + 93.39*fb + 28.52,nfr,1);       % convert to power per erb
0356         dblab='Power/Erb-rate';
0357 end
0358 %
0359 % Now map onto the desired frequency scale
0360 %
0361 if any(mode=='n')
0362     fbopt=['cushn' mdsw(1)];
0363 else
0364     fbopt=['cush' mdsw(1)];
0365 end
0366 b=b*filtbankm(nfrq,fftlen,fftfs,fx(1),fx(end),fbopt)';
0367 
0368 if ~nargout || any(mode=='g') ||  any(lower(mode)=='d')
0369     if numel(db)<2          % find clipping limits
0370         plim=max(b(:))*[0.1^(0.1*db) 1];
0371     else
0372         plim=10.^(0.1*db(1:2));
0373     end
0374     if plim(2)<=0
0375         plim(2)=1;
0376     end
0377     if plim(1)<=0 || plim(1)==plim(2)
0378         plim(1)=0.1*plim(2);
0379     end
0380     if ~nargout || any(mode=='g')
0381         bd=10*log10(max(b,max(b(:)*1e-30)));  % save an unclipped log version for plotting
0382     end
0383     if any(mode=='D')
0384         b=min(max(b,plim(1)),plim(2)); % clip the output
0385     end
0386     if any(mode=='d')
0387         b=10*log10(b);    % output the dB version
0388     end
0389 end
0390 % now plot things
0391 if ~nargout || any(mode=='g')
0392     cla;  % clear current axis
0393     imagesc(t,fx,bd');
0394     axis('xy');
0395     set(gca,'tickdir','out','clim',10*log10(plim));
0396     if any(mode=='j')
0397         colormap('jet');
0398         map=colormap;
0399     elseif any(mode=='J')
0400         map=tcmap;
0401     else
0402         map = repmat((0:63)'/63,1,3);
0403     end
0404     if any(mode=='i')               % 'i' option = invert the colourmap
0405         map=map(64:-1:1,:);
0406     end
0407     colormap(map);
0408     if any(mode=='c')                % 'c' option = show a colourbar
0409         colorbar;
0410         v_cblabel([dblab ' (dB)']);
0411     end
0412     %
0413     % Now check if annotations or a waveform are required
0414     %
0415     dotaw=[((any(mode=='t') && size(ann,2)>1) || size(ann,2)==1) size(ann,2)>1 (any(mode=='w') && ns2==1)];
0416     ylim=get(gca,'ylim');
0417     if  any(dotaw)
0418         yrange = ylim(2)-ylim(1);
0419         zlim=ylim;
0420         toptaw=cumsum([0 dotaw.*[0.05 0.05 0.1]]*yrange)+ylim(2);
0421         zlim(2)=toptaw(4);
0422         set(gca,'ylim',zlim,'color',map(1,:));
0423         if dotaw(3)        % Plot the waveform
0424             six=min(max(floor((get(gca,'xlim')-fs(2))*fs(1))+[1 2],1),ns1);
0425             smax=max(s(six(1):six(2)));
0426             smin=min(s(six(1):six(2)));
0427             if smax==smin
0428                 smax=smax+1;
0429                 smin=smin-1;
0430             end
0431             srange=smax-smin;
0432             hold on
0433             plot(fs(2)+(six(1)-1:six(2)-1)/fs(1),(s(six(1):six(2))-smin)/srange*0.9*(toptaw(4)-toptaw(3))+toptaw(3),'color',map(48,:))
0434             hold off
0435         end
0436         if dotaw(1) || dotaw(2)
0437             tmk=cell2mat(ann(:,1));
0438             tmksel=tmk(:,1)<=t(end) & tmk(:,end)>=t(1);
0439             yix=1+[tmk(tmksel,1)<t(1) ones(sum(tmksel),2) tmk(tmksel,end)>t(end)]';
0440             tmk(:,1)=max(tmk(:,1),t(1));  % clip to axis limits
0441             tmk(:,end)=min(tmk(:,end),t(end));
0442         end
0443         if dotaw(1) && any(tmksel)  % draw time markers
0444             ymk=toptaw(1:2)*[0.8 0.4;0.2 0.6];
0445             switch size(tmk,2)
0446                 case 0
0447                 case 1      % isolated marks
0448                     hold on
0449                     plot([tmk(tmksel) tmk(tmksel)]',repmat(ymk',1,sum(tmksel)),'color',map(48,:));
0450                     hold off
0451                 otherwise % draw durations
0452 
0453                     hold on
0454                     plot(tmk(tmksel,[1 1 2 2])',ymk(yix),'color',map(48,:));
0455                     hold off
0456             end
0457         end
0458         if dotaw(2) && any(tmksel) % print annotations
0459             if any(mode=='a')
0460                 horal='center';
0461                 tmk=(tmk(:,1)+tmk(:,end))*0.5;
0462             else
0463                 horal='left';
0464                 tmk=tmk(:,1);
0465             end
0466             if size(ann,2)>2
0467                 font='Arial';
0468                 for i=1:size(ann,1)
0469                     if tmksel(i)
0470                         if ~isempty(ann{i,3})
0471                             font = ann{i,3};
0472                         end
0473                         text(tmk(i),toptaw(2),ann{i,2},'color',map(48,:),'fontname',font,'VerticalAlignment','baseline','HorizontalAlignment',horal);
0474                     end
0475                 end
0476             else
0477                 for i=1:size(ann,1)
0478                     if tmksel(i)
0479                         text(tmk(i),toptaw(2),ann{i,2},'color',map(48,:),'VerticalAlignment','baseline','HorizontalAlignment',horal);
0480                     end
0481                 end
0482             end
0483         end
0484     end
0485     xlabel(['Time (' v_xticksi 's)']);
0486     if any(mode=='f') && ~strcmp(frlab,'Hz')
0487         ylabel([frlabf '-scaled frequency (Hz)']);
0488         ytickhz(frq2y,y2frq);
0489     else
0490         ylabel(['Frequency (' v_yticksi frlab ')']);
0491     end
0492     ytick=get(gca,'YTick');
0493     ytickl=get(gca,'YTickLabel');
0494     msk=ytick<=ylim(2);
0495     if any(~msk)
0496         set(gca,'YTick',ytick(msk),'YTickLabel',ytickl(msk));
0497     end
0498 end
0499 
0500 function ytickhz(frq2y,y2frq)
0501 % label non linear y frequency axis
0502 %
0503 % Bugs/Suggestions:
0504 % * Add a penalty for large numbers (e.g. 94 is less "round" than 11)
0505 % * possibly add subticks at 1:2:5 if boundaries are 1 and 10
0506 % * could treat subtick allocation specially if bounding lables are both powers of 10
0507 %   and work in log spacing rather than spacing directly
0508 
0509 % algorithm constants
0510 
0511 seps=[0.4 1 3 6]; % spacings: (a) min subtick, (b) min tick, (c) min good tick, (d) max good tick
0512 ww=[0.5 0.6 0.8 0.1 0.3 0.3 0.2];  % weight for (a) last digit=5, (b) power of 10, (c) power of 1000, (d) equal spacing, (e) 1:2:5 labels (f) <seps(3) (g) >seps(4)
0513 nbest=10; % number of possibilities to track
0514 
0515 prefix={'y','z','a','f','p','n','u','m','','k','M','G','T','P','E','Z','Y'};
0516 
0517 ah=gca;
0518 getgca=get(ah);  % Get original axis properties
0519 set(ah,'Units','points','FontUnits','points');
0520 getgcac=get(ah);  % Get axis properties in points units
0521 set(ah,'Units',getgca.Units,'FontUnits',getgca.FontUnits); % return to original values
0522 ylim=getgca.YLim;
0523 yrange=ylim*[-1;1];
0524 chsz= yrange*getgcac.FontSize/getgcac.Position(4); % char height in Y-units
0525 % divide the y-axis up into bins containing at most one label each
0526 maxl=ceil(2*yrange/chsz);  % max number of labels
0527 
0528 % candidate array [cand(:,[1 2])/1000 cand(:,5) cand(:,6)/1000 cand(:,[7 8])]
0529 % 1,2=y limits, 3,4=log limits, 5=Hz, 6=cost, 7=mantissa, 8=exponent, 9=sig digits, 10=y-position
0530 cand=zeros(maxl+2,10);
0531 yinc=(yrange+chsz*0.0002)/maxl;  % bin spacing (allowing for a tiny bit to ensure the ends are included)
0532 cand(2:end-1,2)=ylim(1)+yinc*(1:maxl)'-chsz*0.0001;
0533 cand(3:end-1,1)=cand(2:end-2,2);
0534 cand(2,1)=cand(2,2)-yinc;
0535 cand(2:end-1,1:2)=y2frq(max(cand(2:end-1,1:2),0));
0536 
0537 % find the "roundest" number in each interval
0538 % first deal with intervals containing zero
0539 cand([1 maxl+2],6)=-1;
0540 cand(2,9)=(cand(2,1)<=0);  % mask out interval contaiing zero
0541 cand(2,6)=-cand(2,9);
0542 msk=cand(:,6)==0;  % find rows without a cost yet
0543 cand(msk,3:4)=log10(cand(msk,1:2));
0544 % find powers of 1000
0545 loglim=ceil(cand(:,3:4)/3);
0546 msk=loglim(:,2)>loglim(:,1);
0547 if any(msk)
0548     xp=loglim(msk,1);
0549     wuns=ones(length(xp),1);
0550     cand(msk,5:9)=[1000.^xp wuns-ww(3) wuns 3*xp wuns];
0551 end
0552 % find powers of 10
0553 loglim=ceil(cand(:,3:4));
0554 msk=~msk & (loglim(:,2)>loglim(:,1));
0555 if any(msk)
0556     xp=loglim(msk,1);
0557     wuns=ones(length(xp),1);
0558     cand(msk,5:9)=[10.^xp wuns-ww(2) wuns xp wuns];
0559 end
0560 % find value with fewest digits
0561 msk=cand(:,6)==0;  % find rows without a cost yet
0562 maxsig=1-floor(log10(10^min(cand(msk,3:4)*[-1;1])-1)); % maximum number of significant figures to consider
0563 pten=10.^(0:maxsig-1);   % row vector of powers of ten
0564 noten=10.^(-floor(cand(msk,3))); % exponent of floating point representation of lower bound
0565 sigdig=sum((ceil(cand(msk,2).*noten*pten)-ceil(cand(msk,1).*noten*pten))==0,2); % number of digits common to the interval bounds
0566 lowman=ceil(cand(msk,1).*noten.*10.^sigdig);
0567 midman=10*floor(lowman/10)+5;
0568 highman=ceil(cand(msk,2).*noten.*10.^sigdig);
0569 mskman=midman>=lowman & midman<highman;   % check if we can include a manitssa ending in 5
0570 lowman(mskman)=midman(mskman);
0571 cand(msk,6:9)=[sigdig+1 lowman floor(cand(msk,3))-sigdig sigdig+1];
0572 cand(msk,5)=cand(msk,7).*10.^cand(msk,8);
0573 cand(msk,6)=cand(msk,6)-(mod(cand(msk,7),10)==5)*ww(1);
0574 cand(2:end-1,10)=frq2y(cand(2:end-1,5));
0575 cand([1 maxl+2],10)=ylim + seps(4)*chsz*[-1 1]; % put imaginary labels at the optimum spacing beyond the axes
0576 % [cand(:,[1 2 5])/1000 cand(:,[6 7 8 9])]
0577 
0578 % Now do n-best DP to find the best sequence
0579 
0580 ratint=[8/5 25/10 0 0 4/3];
0581 costs=Inf(nbest,maxl+2); % cumulative path costs
0582 costs(1,1)=0; % starting node only has one option
0583 prev=ones(nbest,maxl+2); % previous label in path
0584 labcnt=zeros(nbest,maxl+2); % number of labels in path
0585 for i=2:maxl+2
0586     ntry=nbest*(i-1); % number of previous options
0587     prevc=reshape(repmat(1:i-1,nbest,1),ntry,1); % previous candidate
0588     prevprev=1+floor((prev(1:ntry)'-1)/nbest); % previous previous candidate
0589     msk=prevprev>1+(maxl+2)*(i==maxl+2); % mask for label triplets
0590     labcnti=labcnt(1:ntry)+1;
0591     disti=(cand(i,10)-cand(prevc,10))/chsz; % distance to previous label in characters
0592     costa=max(seps(3)-disti,0)*ww(6)+max(disti-seps(4),0)*ww(7);
0593     incri=(cand(i,5)-cand(prevc,5)); % label increment
0594     incrj=(cand(i,5)-cand(prevprev,5)); % double label increment
0595     if any(msk)
0596         costa(msk)=costa(msk)- ww(4)*(abs(incrj(msk)-2*incri(msk))<0.01*incri(msk));
0597         if cand(i,7)==1 || cand(i,7)==2 || cand(i,7)==5 % look for labels 1:2:5
0598             costa(msk)=costa(msk)- ww(5)*(abs(incrj(msk)-ratint(cand(i,7))*incri(msk))<0.01*incri(msk));
0599         end
0600     end
0601     costa(disti<seps(2))=Inf;
0602     costi=(costs(1:ntry).*max(labcnt(1:ntry),1)+costa'+cand(i,6))./labcnti;
0603     [sc,isc]=sort(costi);
0604     isc=isc(1:nbest);
0605     costs(:,i)=sc(1:nbest)';
0606     prev(:,i)=isc';
0607     labcnt(:,i)=labcnti(isc)';
0608 end
0609 
0610 % now traceback the best sequence
0611 
0612 % fprintf('Traceback\n\n');
0613 ichoose=0;
0614 labchoose=[];
0615 for i=1:nbest
0616     if labcnt(i,maxl+2)>1 && costs(i,maxl+2)<Inf
0617         lablist=zeros(labcnt(i,maxl+2)-1,1);
0618         k=prev(i,maxl+2);
0619         for j=labcnt(i,maxl+2)-1:-1:1
0620             lablist(j)=1+floor((k-1)/nbest);
0621             k=prev(k);
0622         end
0623         %         fprintf('Cost=%8.2f :',costs(i,maxl+2));
0624         %         fprintf(' %g',cand(lablist,5))
0625         %         fprintf('\n');
0626         if ~ichoose || labcnt(ichoose,maxl+2)==1
0627             ichoose=i;
0628             labchoose=lablist;
0629         end
0630     end
0631 end
0632 
0633 % now create the labels
0634 
0635 ntick=length(labchoose);
0636 % sort out the subticks
0637 subpos=[];
0638 if ntick>=2
0639     for i=1:ntick-1
0640         clj=cand(labchoose(i:i+1),:);
0641         sprec=min(clj(1,8)+100*(clj(1,7)==0),clj(2,8)); % subtick precision
0642         spos=(clj(1,7)*10^(clj(1,8)-sprec):clj(2,7)*10^(clj(2,8)-sprec))*10^sprec;
0643         nsub=length(spos);
0644         if nsub==2
0645             spos=spos*[1 0.5 0;0 0.5 1];
0646             nsub=3;
0647         end
0648         if nsub>=3
0649             yspos=frq2y(spos);
0650             for kk=1:3 % try various subdivisions: every 1, 2 or 5
0651                 k=kk+2*(kk==3);  % 1, 2 and 5
0652                 if 2*k<=nsub-1 && ~mod(nsub-1,k)  % must divide exactly into nsub
0653                     if all((yspos(1+k:k:nsub)-yspos(1:k:nsub-k))>=(seps(1)*chsz)) % check they all fit in
0654                         subpos=[subpos yspos(1+k:k:nsub-k)];
0655                         if i==1
0656                             spos=(ceil(cand(2,1)/10^sprec):clj(1,7)*10^(clj(1,8)-sprec))*10^sprec;
0657                             nsub=length(spos);
0658                             yspos=frq2y(spos);
0659                             if nsub>=k+1 && all((yspos(nsub:-k:1+k)-yspos(nsub-k:-k:1))>=(seps(1)*chsz))
0660                                 subpos=[subpos yspos(nsub-k:-k:1)];
0661                             end
0662                         elseif i==ntick-1
0663                             spos=(clj(2,7)*10^(clj(2,8)-sprec):floor(cand(end-1,2)/10^sprec))*10^sprec;
0664                             nsub=length(spos);
0665                             yspos=frq2y(spos);
0666                             if nsub>=k+1 && all((yspos(1+k:k:nsub)-yspos(1:k:nsub-k))>=(seps(1)*chsz))
0667                                 subpos=[subpos yspos(1+k:k:nsub)];
0668                             end
0669                         end
0670                         break;
0671                     end
0672                 end
0673             end
0674         end
0675     end
0676 end
0677 nsub=length(subpos);
0678 tickpos=[cand(labchoose,10); subpos'];
0679 ticklab=cell(ntick+nsub,1);
0680 sipref=min(max(floor((sum(cand(labchoose,8:9),2)-1)/3),-8),8);
0681 nzadd=cand(labchoose,8)-3*sipref;  % trailing zeros to add
0682 digzer=cand(labchoose,7).*10.^max(nzadd,0); % label digits including trailing zeros
0683 ndleft=cand(labchoose,9)+nzadd; % digits to the left of the decimal point
0684 for i=1:ntick
0685     tickint=num2str(digzer(i));
0686     if nzadd(i)<0
0687         tickint=[tickint(1:ndleft(i)) '.' tickint(1+ndleft(i):end)];
0688     end
0689     ticklab{i} = sprintf('%s%s',tickint,prefix{sipref(i)+9});
0690 end
0691 for i=ntick+1:ntick+nsub
0692     ticklab{i}='';
0693 end
0694 [tickpos,ix]=sort(tickpos);
0695 ticklab=ticklab(ix);
0696 
0697 set(ah,'YTick',tickpos','YTickLabel',ticklab);
0698

Generated by m2html © 2003