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 % * 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: v_spgrambw.m 10865 2018-09-21 17:22:45Z 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,v_frq2mel(fnyq),nfrq); % DC to Nyquist 0238 case 'b' 0239 fx=linspace(0,v_frq2bark(fnyq),nfrq); % DC to Nyquist 0240 case 'e' 0241 fx=linspace(0,v_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=v_frq2mel(fmax); % DC to Nyquist 0252 case 'b' 0253 fmaxu=v_frq2bark(fmax); % DC to Nyquist 0254 case 'e' 0255 fmaxu=v_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=v_mel2frq(fx); 0285 frlab='Mel'; 0286 frlabf='Mel'; 0287 frq2y=@v_frq2mel; 0288 y2frq=@v_mel2frq; 0289 case 'b' 0290 f=v_bark2frq(fx); 0291 frlab='Bark'; 0292 frlabf='Bark'; 0293 frq2y=@v_frq2bark; 0294 y2frq=@v_bark2frq; 0295 case 'e' 0296 f=v_erb2frq(fx); 0297 frlab='Erb-rate'; 0298 frlabf='Erb'; 0299 frq2y=@v_frq2erb; 0300 y2frq=@v_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]=v_enframe(s(ix1:ix2),win,ninc); 0323 t=fs(2)+(t+ix1-2)/fs(1); % time axis 0324 b=v_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 v_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 v_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 (' v_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 (' v_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','u','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