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