Home > voicebox > spgrambw.m

spgrambw

PURPOSE ^

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

SYNOPSIS ^

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

DESCRIPTION ^

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

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

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

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

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

Generated on Fri 22-Sep-2017 19:37:38 by m2html © 2003