Description of spgrambw

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

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
