# 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 word transcription

             (3) [s,fs,wrd,phn]=v_readsph(filename);           % read a TIMIT file
                 v_spgrambw(s,fs,'pJcwaAtT',[],[],[],[],{wrd phn}); % plot spectrogram with word and phone transcriptions

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

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

             (6) 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 (Helvetica) will be used (To display: g=get(gca);g.Children(1).FontName).
                         For multiple annotation lines, ANN should be a cell  array whose elements are cell arrays as described above.

     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','A' = centre-align annotations rather than left-aligning them.
                      'a' applies to the lowest annotation, 'A' to all the others.
            't','T' = add time markers with annotations.
                      't' applies to the lowest annotation, 'T' to all the others.

     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:

- v_bark2frq V_BARK2FRQ Convert the BARK frequency scale to Hertz FRQ=(BARK)
- v_cblabel V_CBLABEL add a label to a colorbar C=(L,H)
- v_enframe V_ENFRAME split signal up into (overlapping) frames: one per row. [F,T]=(X,WIN,HOP)
- v_erb2frq V_ERB2FRQ Convert ERB frequency scale to Hertz FRQ=(ERB)
- v_frq2bark V_FRQ2BARK Convert Hertz to BARK frequency scale BARK=(FRQ)
- v_frq2erb V_FRQ2ERB Convert Hertz to ERB frequency scale ERB=(FRQ)
- v_frq2mel V_FRQ2ERB Convert Hertz to Mel frequency scale MEL=(FRQ)
- v_mel2frq V_MEL2FRQ Convert Mel frequency scale to Hertz FRQ=(MEL)
- v_rfft V_RFFT Calculate the DFT of real data Y=(X,N,D)
- v_xticksi V_XTIXKSI labels the x-axis of a plot using SI multipliers S=(AH)
- v_yticksi V_YTIXKSI labels the y-axis of a plot using SI multipliers S=(AH)

This function is called by:

- v_activlevg V_ACTIVLEVG Measure active speech level robustly [LEV,AF,FSO]=(sp,FS,MODE)
- v_fxpefac V_FXPEFAC PEFAC pitch tracker [FX,TT,PV,FV]=(S,FS,TINC,M,PP)
- v_readwav V_READWAV Read a .WAV format sound file [Y,FS,WMODE,FIDX]=(FILENAME,MODE,NMAX,NSKIP)
- v_stdspectrum V_STDSPECTRUM Generate standard acoustic/speech spectra in s- or z-domain [B,A,SI,SN]=(S,M,F,N,ZI,BS,AS)

## SUBFUNCTIONS

- function ytickhz(frq2y,y2frq)

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