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) [s,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) [s,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 %
0142 % first decode the input arguments
0143 %
0144 if size(s,1)==1
0145     s=s(:);   % force to be a column vector (unless it is a matrix)
0146 end
0147 [ns1,ns2]=size(s);
0148 ap=zeros(1,6);
0149 j=2;
0150 if numel(fs)<2
0151     fs(2)=1/fs(1);  % first sample or frame is at time 1/fs
0152 end
0153 for i=1:length(varargin)
0154     if ischar(varargin{i})
0155         ap(1)=i;
0156     else
0157         ap(j)=i;
0158         j=j+1;
0159     end
0160 end
0161 if ap(1) && ~isempty(varargin{ap(1)})
0162     mode=varargin{ap(1)};
0163 else
0164     mode='';  % default mode
0165 end
0166 if ap(2) && ~isempty(varargin{ap(2)})
0167     bw=varargin{ap(2)};
0168 else
0169     bw=200;
0170 end
0171 if ap(3) && ~isempty(varargin{ap(3)})
0172     fmax=varargin{ap(3)};
0173 else
0174     fmax=[];
0175 end
0176 if ap(4) && ~isempty(varargin{ap(4)})
0177     db=varargin{ap(4)};
0178 else
0179     db=40;
0180 end
0181 if ap(5) && ~isempty(varargin{ap(5)})
0182     tinc=varargin{ap(5)};
0183 else
0184     tinc=0;
0185 end
0186 switch numel(tinc)
0187     case 1
0188         tinc=[tinc -Inf Inf];
0189     case 2
0190         tinc=[0 tinc];
0191     otherwise
0192         tinc=tinc([2 1 3]);
0193 end
0194 if tinc(1)<=0
0195     tinc(1)=1.81/(4*bw); % default frame increment
0196 end
0197 if ap(6)
0198     ann=varargin{ap(6)};
0199 else
0200     ann=[];
0201 end
0202 
0203 % now sort out the mode flags
0204 
0205 mdsw='  ';           % [yscale preemph]
0206 for i=1:length(mode)
0207     switch mode(i)
0208         case {'l','m','b','e'}
0209             mdsw(1)=mode(i);
0210         case {'p','P'}
0211             mdsw(2)=mode(i);
0212     end
0213 end
0214 if mdsw(2)=='P'
0215     mdsw(2)=mdsw(1);        % preemphasis is scaling dependent
0216 end
0217 %
0218 % sort out the frequency axis
0219 %
0220 flmin=30;                   % min frequency for 'l' option
0221 nfrq=257;                   % default number of frequency bins
0222 if ns2==1
0223     fnyq=fs(1)/2;           % default upper frequency limit is fs/2
0224 else                        % input is a power spectrum
0225     if numel(fs)<3
0226         fs(3)=fs(1)*0.25;   % default increment is 0.25 times frame increment
0227     end
0228     if numel(fs)<4
0229         fs(4)=0;            % first freq bin is DC by default
0230     end
0231     fnyq=fs(4)+(ns2-1)*fs(3);  % default upper frequency limit is highest supplied frequency
0232 end
0233 
0234 if ~numel(fmax)             % no explicit frequency range
0235     switch mdsw(1)
0236         case 'l'
0237             fx=linspace(log10(flmin),log10(fnyq),nfrq);   % 20  Hz to Nyquist
0238         case 'm'
0239             fx=linspace(0,v_frq2mel(fnyq),nfrq);   % DC to Nyquist
0240         case 'b'
0241             fx=linspace(0,v_frq2bark(fnyq),nfrq);   % DC to Nyquist
0242         case 'e'
0243             fx=linspace(0,v_frq2erb(fnyq),nfrq);   % DC to Nyquist
0244         otherwise   % linear Hz scale
0245             fx=(0:nfrq-1)*fnyq/(nfrq-1);
0246     end
0247 else
0248     if any(mode=='h')
0249         switch mdsw(1)
0250             case 'l'
0251                 fmaxu=log10(fmax);   % 20  Hz to Nyquist
0252             case 'm'
0253                 fmaxu=v_frq2mel(fmax);   % DC to Nyquist
0254             case 'b'
0255                 fmaxu=v_frq2bark(fmax);   % DC to Nyquist
0256             case 'e'
0257                 fmaxu=v_frq2erb(fmax);   % DC to Nyquist
0258             otherwise
0259                 fmaxu=fmax;  % linear Hz scale
0260         end
0261     else
0262         fmaxu=fmax;                 % already in the correct units
0263     end
0264     if numel(fmax)<2   % only max value specified
0265         if mdsw(1)=='l'
0266             fx=linspace(log10(flmin),fmaxu,nfrq);   % 20  Hz to fmax
0267         else
0268             fx=linspace(0,fmaxu,nfrq);   % DC to fmax
0269         end
0270     elseif numel(fmax)<3 % min and max values specified
0271         fx=linspace(fmaxu(1),fmaxu(2),nfrq);   % fmin to fmax
0272     else
0273         fmaxu(2)=fmax(2)*(fmaxu(3)-fmaxu(1))/(fmax(3)-fmax(1)); % scale the step size appropriately
0274         fx=fmaxu(1):fmaxu(2):fmaxu(3);   % fmin to fmax in steps of finc
0275         nfrq=length(fx);
0276     end
0277 end
0278 switch mdsw(1)          % convert the frequency range to Hz
0279     case 'l'
0280         f=10.^fx;
0281         frlab='log_{10}Hz';
0282         frlabf='log';
0283         frq2y=@log10;
0284         y2frq=@(x) 10.^x;
0285     case 'm'
0286         f=v_mel2frq(fx);
0287         frlab='Mel';
0288         frlabf='Mel';
0289         frq2y=@v_frq2mel;
0290         y2frq=@v_mel2frq;
0291     case 'b'
0292         f=v_bark2frq(fx);
0293         frlab='Bark';
0294         frlabf='Bark';
0295         frq2y=@v_frq2bark;
0296         y2frq=@v_bark2frq;
0297     case 'e'
0298         f=v_erb2frq(fx);
0299         frlab='Erb-rate';
0300         frlabf='Erb';
0301         frq2y=@v_frq2erb;
0302         y2frq=@v_erb2frq;
0303     otherwise
0304         f=fx;
0305         frlab='Hz';
0306         frq2y=@(x) x;
0307         y2frq=@(x) x;
0308 end
0309 if ~any(mode=='H')
0310     f=fx;               % give output frequencies in native units instead of Hz unless 'H' is specified
0311 end
0312 %
0313 % now calculate the spectrogram
0314 %
0315 if ns2==1   % input is a speech signal vector
0316     winlen = fix(1.81*fs(1)/bw);   % window length
0317     win=0.54+0.46*cos((1-winlen:2:winlen)*pi/winlen);  % Hamming window
0318     ninc=max(round(tinc(1)*fs(1)),1);                 % window increment in samples
0319     %  we need to take account of minimum freq increment + make it exact if possible
0320     fftlen=pow2(nextpow2(4*winlen));        % enough oversampling to get good interpolation
0321     win=win/sqrt(sum(win.^2));              % ensure window squared sums to unity
0322     ix1=max(round((tinc(2)-fs(2))*fs(1)-(winlen-3)/2),1); % first sample required
0323     ix2=min(ceil((tinc(3)-fs(2))*fs(1)+(winlen+1)/2),ns1); % last sample required
0324     [sf,t]=v_enframe(s(ix1:ix2),win,ninc);
0325     t=fs(2)+(t+ix1-2)/fs(1);                         % time axis
0326     b=v_rfft(sf,fftlen,2);
0327     b=b.*conj(b)*2/fs(1);          % Power per Hz
0328     b(:,1)=b(:,1)*0.5;   % correct for no negative zero frequency to double the power
0329     b(:,end)=b(:,end)*0.5;   % correct for no negative nyquist frequency to double the power
0330     fb=(0:fftlen/2)*fs(1)/fftlen; % fft bin frequencies
0331     fftfs=fs(1);
0332 else
0333     b=s;
0334     t=fs(2)+(0:ns1-1)/fs(1);  % frame times
0335     fb=fs(4)+(0:ns2-1)*fs(3);
0336     fftlen=[ns2 fs(3) fs(4)]; % for v_filtbankm: ns2=# input freq bins, freq increment (Hz), first bin freq (Hz)
0337     fftfs=0;
0338     %     fftlen=2*(ns2-1);  % assume an even length fft
0339     %     fftfs=fftlen*fs(3);
0340 end
0341 nfr=numel(t);                   % number of frames
0342 dblab='Power/Hz';
0343 switch mdsw(2)
0344     case {'p','l'}
0345         b=b.*repmat(fb*log(10),nfr,1);       % convert to power per decade
0346         dblab='Power/Decade';
0347     case 'm'
0348         b=b.*repmat((700+fb)*log(1+1000/700)/1000,nfr,1);       % convert to power per mel
0349         dblab='Power/Mel';
0350     case 'b'
0351         b=b.*repmat((1960+fb).^2/52547.6,nfr,1);       % convert to power per bark
0352         dblab='Power/Bark';
0353     case 'e'
0354         b=b.*repmat(6.23*fb.^2 + 93.39*fb + 28.52,nfr,1);       % convert to power per erb
0355         dblab='Power/Erb-rate';
0356 end
0357 %
0358 % Now map onto the desired frequency scale
0359 %
0360 if any(mode=='n')
0361     fbopt=['cushn' mdsw(1)];
0362 else
0363     fbopt=['cush' mdsw(1)];
0364 end
0365 b=b*filtbankm(nfrq,fftlen,fftfs,fx(1),fx(end),fbopt)';
0366 
0367 if ~nargout || any(mode=='g') ||  any(lower(mode)=='d')
0368     if numel(db)<2          % find clipping limits
0369         plim=max(b(:))*[0.1^(0.1*db) 1];
0370     else
0371         plim=10.^(0.1*db(1:2));
0372     end
0373     if plim(2)<=0
0374         plim(2)=1;
0375     end
0376     if plim(1)<=0 || plim(1)==plim(2)
0377         plim(1)=0.1*plim(2);
0378     end
0379     if ~nargout || any(mode=='g')
0380         bd=10*log10(max(b,max(b(:)*1e-30)));  % save an unclipped log version for plotting
0381     end
0382     if any(mode=='D')
0383         b=min(max(b,plim(1)),plim(2)); % clip the output
0384     end
0385     if any(mode=='d')
0386         b=10*log10(b);    % output the dB version
0387     end
0388 end
0389 % now plot things
0390 if ~nargout || any(mode=='g')
0391     cla;  % clear current axis
0392     imagesc(t,fx,bd');
0393     axis('xy');
0394     set(gca,'tickdir','out','clim',10*log10(plim));
0395     if any(mode=='j')
0396         colormap('jet');
0397         map=colormap;
0398     elseif any(mode=='J')
0399         map=tcmap;
0400     else
0401         map = repmat((0:63)'/63,1,3);
0402     end
0403     if any(mode=='i')               % 'i' option = invert the colourmap
0404         map=map(64:-1:1,:);
0405     end
0406     colormap(map);
0407     if any(mode=='c')                % 'c' option = show a colourbar
0408         colorbar;
0409         v_cblabel([dblab ' (dB)']);
0410     end
0411     %
0412     % Now check if annotations or a waveform are required
0413     %
0414     dotaw=[((any(mode=='t') && size(ann,2)>1) || size(ann,2)==1) size(ann,2)>1 (any(mode=='w') && ns2==1)];
0415     ylim=get(gca,'ylim');
0416     if  any(dotaw)
0417         yrange = ylim(2)-ylim(1);
0418         zlim=ylim;
0419         toptaw=cumsum([0 dotaw.*[0.05 0.05 0.1]]*yrange)+ylim(2);
0420         zlim(2)=toptaw(4);
0421         set(gca,'ylim',zlim,'color',map(1,:));
0422         if dotaw(3)        % Plot the waveform
0423             six=min(max(floor((get(gca,'xlim')-fs(2))*fs(1))+[1 2],1),ns1);
0424             smax=max(s(six(1):six(2)));
0425             smin=min(s(six(1):six(2)));
0426             if smax==smin
0427                 smax=smax+1;
0428                 smin=smin-1;
0429             end
0430             srange=smax-smin;
0431             hold on
0432             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,:))
0433             hold off
0434         end
0435         if dotaw(1) || dotaw(2)
0436             tmk=cell2mat(ann(:,1));
0437             tmksel=tmk(:,1)<=t(end) & tmk(:,end)>=t(1);
0438             yix=1+[tmk(tmksel,1)<t(1) ones(sum(tmksel),2) tmk(tmksel,end)>t(end)]';
0439             tmk(:,1)=max(tmk(:,1),t(1));  % clip to axis limits
0440             tmk(:,end)=min(tmk(:,end),t(end));
0441         end
0442         if dotaw(1) && any(tmksel)  % draw time markers
0443             ymk=toptaw(1:2)*[0.8 0.4;0.2 0.6];
0444             switch size(tmk,2)
0445                 case 0
0446                 case 1      % isolated marks
0447                     hold on
0448                     plot([tmk(tmksel) tmk(tmksel)]',repmat(ymk',1,sum(tmksel)),'color',map(48,:));
0449                     hold off
0450                 otherwise % draw durations
0451 
0452                     hold on
0453                     plot(tmk(tmksel,[1 1 2 2])',ymk(yix),'color',map(48,:));
0454                     hold off
0455             end
0456         end
0457         if dotaw(2) && any(tmksel) % print annotations
0458             if any(mode=='a')
0459                 horal='center';
0460                 tmk=(tmk(:,1)+tmk(:,end))*0.5;
0461             else
0462                 horal='left';
0463                 tmk=tmk(:,1);
0464             end
0465             if size(ann,2)>2
0466                 font='Arial';
0467                 for i=1:size(ann,1)
0468                     if tmksel(i)
0469                         if ~isempty(ann{i,3})
0470                             font = ann{i,3};
0471                         end
0472                         text(tmk(i),toptaw(2),ann{i,2},'color',map(48,:),'fontname',font,'VerticalAlignment','baseline','HorizontalAlignment',horal);
0473                     end
0474                 end
0475             else
0476                 for i=1:size(ann,1)
0477                     if tmksel(i)
0478                         text(tmk(i),toptaw(2),ann{i,2},'color',map(48,:),'VerticalAlignment','baseline','HorizontalAlignment',horal);
0479                     end
0480                 end
0481             end
0482         end
0483     end
0484     xlabel(['Time (' v_xticksi 's)']);
0485     if any(mode=='f') && ~strcmp(frlab,'Hz')
0486         ylabel([frlabf '-scaled frequency (Hz)']);
0487         ytickhz(frq2y,y2frq);
0488     else
0489         ylabel(['Frequency (' v_yticksi frlab ')']);
0490     end
0491     ytick=get(gca,'YTick');
0492     ytickl=get(gca,'YTickLabel');
0493     msk=ytick<=ylim(2);
0494     if any(~msk)
0495         set(gca,'YTick',ytick(msk),'YTickLabel',ytickl(msk));
0496     end
0497 end
0498 
0499 function ytickhz(frq2y,y2frq)
0500 % label non linear y frequency axis
0501 %
0502 % Bugs/Suggestions:
0503 % * Add a penalty for large numbers (e.g. 94 is less "round" than 11)
0504 % * possibly add subticks at 1:2:5 if boundaries are 1 and 10
0505 % * could treat subtick allocation specially if bounding lables are both powers of 10
0506 %   and work in log spacing rather than spacing directly
0507 
0508 % algorithm constants
0509 
0510 seps=[0.4 1 3 6]; % spacings: (a) min subtick, (b) min tick, (c) min good tick, (d) max good tick
0511 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)
0512 nbest=10; % number of possibilities to track
0513 
0514 prefix={'y','z','a','f','p','n','u','m','','k','M','G','T','P','E','Z','Y'};
0515 
0516 ah=gca;
0517 getgca=get(ah);  % Get original axis properties
0518 set(ah,'Units','points','FontUnits','points');
0519 getgcac=get(ah);  % Get axis properties in points units
0520 set(ah,'Units',getgca.Units,'FontUnits',getgca.FontUnits); % return to original values
0521 ylim=getgca.YLim;
0522 yrange=ylim*[-1;1];
0523 chsz= yrange*getgcac.FontSize/getgcac.Position(4); % char height in Y-units
0524 % divide the y-axis up into bins containing at most one label each
0525 maxl=ceil(2*yrange/chsz);  % max number of labels
0526 
0527 % candidate array [cand(:,[1 2])/1000 cand(:,5) cand(:,6)/1000 cand(:,[7 8])]
0528 % 1,2=y limits, 3,4=log limits, 5=Hz, 6=cost, 7=mantissa, 8=exponent, 9=sig digits, 10=y-position
0529 cand=zeros(maxl+2,10);
0530 yinc=(yrange+chsz*0.0002)/maxl;  % bin spacing (allowing for a tiny bit to ensure the ends are included)
0531 cand(2:end-1,2)=ylim(1)+yinc*(1:maxl)'-chsz*0.0001;
0532 cand(3:end-1,1)=cand(2:end-2,2);
0533 cand(2,1)=cand(2,2)-yinc;
0534 cand(2:end-1,1:2)=y2frq(max(cand(2:end-1,1:2),0));
0535 
0536 % find the "roundest" number in each interval
0537 % first deal with intervals containing zero
0538 cand([1 maxl+2],6)=-1;
0539 cand(2,9)=(cand(2,1)<=0);  % mask out interval contaiing zero
0540 cand(2,6)=-cand(2,9);
0541 msk=cand(:,6)==0;  % find rows without a cost yet
0542 cand(msk,3:4)=log10(cand(msk,1:2));
0543 % find powers of 1000
0544 loglim=ceil(cand(:,3:4)/3);
0545 msk=loglim(:,2)>loglim(:,1);
0546 if any(msk)
0547     xp=loglim(msk,1);
0548     wuns=ones(length(xp),1);
0549     cand(msk,5:9)=[1000.^xp wuns-ww(3) wuns 3*xp wuns];
0550 end
0551 % find powers of 10
0552 loglim=ceil(cand(:,3:4));
0553 msk=~msk & (loglim(:,2)>loglim(:,1));
0554 if any(msk)
0555     xp=loglim(msk,1);
0556     wuns=ones(length(xp),1);
0557     cand(msk,5:9)=[10.^xp wuns-ww(2) wuns xp wuns];
0558 end
0559 % find value with fewest digits
0560 msk=cand(:,6)==0;  % find rows without a cost yet
0561 maxsig=1-floor(log10(10^min(cand(msk,3:4)*[-1;1])-1)); % maximum number of significant figures to consider
0562 pten=10.^(0:maxsig-1);   % row vector of powers of ten
0563 noten=10.^(-floor(cand(msk,3))); % exponent of floating point representation of lower bound
0564 sigdig=sum((ceil(cand(msk,2).*noten*pten)-ceil(cand(msk,1).*noten*pten))==0,2); % number of digits common to the interval bounds
0565 lowman=ceil(cand(msk,1).*noten.*10.^sigdig);
0566 midman=10*floor(lowman/10)+5;
0567 highman=ceil(cand(msk,2).*noten.*10.^sigdig);
0568 mskman=midman>=lowman & midman<highman;   % check if we can include a manitssa ending in 5
0569 lowman(mskman)=midman(mskman);
0570 cand(msk,6:9)=[sigdig+1 lowman floor(cand(msk,3))-sigdig sigdig+1];
0571 cand(msk,5)=cand(msk,7).*10.^cand(msk,8);
0572 cand(msk,6)=cand(msk,6)-(mod(cand(msk,7),10)==5)*ww(1);
0573 cand(2:end-1,10)=frq2y(cand(2:end-1,5));
0574 cand([1 maxl+2],10)=ylim + seps(4)*chsz*[-1 1]; % put imaginary labels at the optimum spacing beyond the axes
0575 % [cand(:,[1 2 5])/1000 cand(:,[6 7 8 9])]
0576 
0577 % Now do n-best DP to find the best sequence
0578 
0579 ratint=[8/5 25/10 0 0 4/3];
0580 costs=Inf(nbest,maxl+2); % cumulative path costs
0581 costs(1,1)=0; % starting node only has one option
0582 prev=ones(nbest,maxl+2); % previous label in path
0583 labcnt=zeros(nbest,maxl+2); % number of labels in path
0584 for i=2:maxl+2
0585     ntry=nbest*(i-1); % number of previous options
0586     prevc=reshape(repmat(1:i-1,nbest,1),ntry,1); % previous candidate
0587     prevprev=1+floor((prev(1:ntry)'-1)/nbest); % previous previous candidate
0588     msk=prevprev>1+(maxl+2)*(i==maxl+2); % mask for label triplets
0589     labcnti=labcnt(1:ntry)+1;
0590     disti=(cand(i,10)-cand(prevc,10))/chsz; % distance to previous label in characters
0591     costa=max(seps(3)-disti,0)*ww(6)+max(disti-seps(4),0)*ww(7);
0592     incri=(cand(i,5)-cand(prevc,5)); % label increment
0593     incrj=(cand(i,5)-cand(prevprev,5)); % double label increment
0594     if any(msk)
0595         costa(msk)=costa(msk)- ww(4)*(abs(incrj(msk)-2*incri(msk))<0.01*incri(msk));
0596         if cand(i,7)==1 || cand(i,7)==2 || cand(i,7)==5 % look for labels 1:2:5
0597             costa(msk)=costa(msk)- ww(5)*(abs(incrj(msk)-ratint(cand(i,7))*incri(msk))<0.01*incri(msk));
0598         end
0599     end
0600     costa(disti<seps(2))=Inf;
0601     costi=(costs(1:ntry).*max(labcnt(1:ntry),1)+costa'+cand(i,6))./labcnti;
0602     [sc,isc]=sort(costi);
0603     isc=isc(1:nbest);
0604     costs(:,i)=sc(1:nbest)';
0605     prev(:,i)=isc';
0606     labcnt(:,i)=labcnti(isc)';
0607 end
0608 
0609 % now traceback the best sequence
0610 
0611 % fprintf('Traceback\n\n');
0612 ichoose=0;
0613 labchoose=[];
0614 for i=1:nbest
0615     if labcnt(i,maxl+2)>1 && costs(i,maxl+2)<Inf
0616         lablist=zeros(labcnt(i,maxl+2)-1,1);
0617         k=prev(i,maxl+2);
0618         for j=labcnt(i,maxl+2)-1:-1:1
0619             lablist(j)=1+floor((k-1)/nbest);
0620             k=prev(k);
0621         end
0622         %         fprintf('Cost=%8.2f :',costs(i,maxl+2));
0623         %         fprintf(' %g',cand(lablist,5))
0624         %         fprintf('\n');
0625         if ~ichoose || labcnt(ichoose,maxl+2)==1
0626             ichoose=i;
0627             labchoose=lablist;
0628         end
0629     end
0630 end
0631 
0632 % now create the labels
0633 
0634 ntick=length(labchoose);
0635 % sort out the subticks
0636 subpos=[];
0637 if ntick>=2
0638     for i=1:ntick-1
0639         clj=cand(labchoose(i:i+1),:);
0640         sprec=min(clj(1,8)+100*(clj(1,7)==0),clj(2,8)); % subtick precision
0641         spos=(clj(1,7)*10^(clj(1,8)-sprec):clj(2,7)*10^(clj(2,8)-sprec))*10^sprec;
0642         nsub=length(spos);
0643         if nsub==2
0644             spos=spos*[1 0.5 0;0 0.5 1];
0645             nsub=3;
0646         end
0647         if nsub>=3
0648             yspos=frq2y(spos);
0649             for kk=1:3 % try various subdivisions: every 1, 2 or 5
0650                 k=kk+2*(kk==3);  % 1, 2 and 5
0651                 if 2*k<=nsub-1 && ~mod(nsub-1,k)  % must divide exactly into nsub
0652                     if all((yspos(1+k:k:nsub)-yspos(1:k:nsub-k))>=(seps(1)*chsz)) % check they all fit in
0653                         subpos=[subpos yspos(1+k:k:nsub-k)];
0654                         if i==1
0655                             spos=(ceil(cand(2,1)/10^sprec):clj(1,7)*10^(clj(1,8)-sprec))*10^sprec;
0656                             nsub=length(spos);
0657                             yspos=frq2y(spos);
0658                             if nsub>=k+1 && all((yspos(nsub:-k:1+k)-yspos(nsub-k:-k:1))>=(seps(1)*chsz))
0659                                 subpos=[subpos yspos(nsub-k:-k:1)];
0660                             end
0661                         elseif i==ntick-1
0662                             spos=(clj(2,7)*10^(clj(2,8)-sprec):floor(cand(end-1,2)/10^sprec))*10^sprec;
0663                             nsub=length(spos);
0664                             yspos=frq2y(spos);
0665                             if nsub>=k+1 && all((yspos(1+k:k:nsub)-yspos(1:k:nsub-k))>=(seps(1)*chsz))
0666                                 subpos=[subpos yspos(1+k:k:nsub)];
0667                             end
0668                         end
0669                         break;
0670                     end
0671                 end
0672             end
0673         end
0674     end
0675 end
0676 nsub=length(subpos);
0677 tickpos=[cand(labchoose,10); subpos'];
0678 ticklab=cell(ntick+nsub,1);
0679 sipref=min(max(floor((sum(cand(labchoose,8:9),2)-1)/3),-8),8);
0680 nzadd=cand(labchoose,8)-3*sipref;  % trailing zeros to add
0681 digzer=cand(labchoose,7).*10.^max(nzadd,0); % label digits including trailing zeros
0682 ndleft=cand(labchoose,9)+nzadd; % digits to the left of the decimal point
0683 for i=1:ntick
0684     tickint=num2str(digzer(i));
0685     if nzadd(i)<0
0686         tickint=[tickint(1:ndleft(i)) '.' tickint(1+ndleft(i):end)];
0687     end
0688     ticklab{i} = sprintf('%s%s',tickint,prefix{sipref(i)+9});
0689 end
0690 for i=ntick+1:ntick+nsub
0691     ticklab{i}='';
0692 end
0693 [tickpos,ix]=sort(tickpos);
0694 ticklab=ticklab(ix);
0695 
0696 set(ah,'YTick',tickpos','YTickLabel',ticklab);
0697

Generated by m2html © 2003