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