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) [s,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) [s,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 % 0142 % first decode the input arguments 0143 % 0144 if size(s,1)==1 0145 s=s(:); % force to be a column vector (unless it is a matrix) 0146 end 0147 [ns1,ns2]=size(s); 0148 ap=zeros(1,6); 0149 j=2; 0150 if numel(fs)<2 0151 fs(2)=1/fs(1); % first sample or frame is at time 1/fs 0152 end 0153 for i=1:length(varargin) 0154 if ischar(varargin{i}) 0155 ap(1)=i; 0156 else 0157 ap(j)=i; 0158 j=j+1; 0159 end 0160 end 0161 if ap(1) && ~isempty(varargin{ap(1)}) 0162 mode=varargin{ap(1)}; 0163 else 0164 mode=''; % default mode 0165 end 0166 if ap(2) && ~isempty(varargin{ap(2)}) 0167 bw=varargin{ap(2)}; 0168 else 0169 bw=200; 0170 end 0171 if ap(3) && ~isempty(varargin{ap(3)}) 0172 fmax=varargin{ap(3)}; 0173 else 0174 fmax=[]; 0175 end 0176 if ap(4) && ~isempty(varargin{ap(4)}) 0177 db=varargin{ap(4)}; 0178 else 0179 db=40; 0180 end 0181 if ap(5) && ~isempty(varargin{ap(5)}) 0182 tinc=varargin{ap(5)}; 0183 else 0184 tinc=0; 0185 end 0186 switch numel(tinc) 0187 case 1 0188 tinc=[tinc -Inf Inf]; 0189 case 2 0190 tinc=[0 tinc]; 0191 otherwise 0192 tinc=tinc([2 1 3]); 0193 end 0194 if tinc(1)<=0 0195 tinc(1)=1.81/(4*bw); % default frame increment 0196 end 0197 if ap(6) 0198 ann=varargin{ap(6)}; 0199 else 0200 ann=[]; 0201 end 0202 0203 % now sort out the mode flags 0204 0205 mdsw=' '; % [yscale preemph] 0206 for i=1:length(mode) 0207 switch mode(i) 0208 case {'l','m','b','e'} 0209 mdsw(1)=mode(i); 0210 case {'p','P'} 0211 mdsw(2)=mode(i); 0212 end 0213 end 0214 if mdsw(2)=='P' 0215 mdsw(2)=mdsw(1); % preemphasis is scaling dependent 0216 end 0217 % 0218 % sort out the frequency axis 0219 % 0220 flmin=30; % min frequency for 'l' option 0221 nfrq=257; % default number of frequency bins 0222 if ns2==1 0223 fnyq=fs(1)/2; % default upper frequency limit is fs/2 0224 else % input is a power spectrum 0225 if numel(fs)<3 0226 fs(3)=fs(1)*0.25; % default increment is 0.25 times frame increment 0227 end 0228 if numel(fs)<4 0229 fs(4)=0; % first freq bin is DC by default 0230 end 0231 fnyq=fs(4)+(ns2-1)*fs(3); % default upper frequency limit is highest supplied frequency 0232 end 0233 0234 if ~numel(fmax) % no explicit frequency range 0235 switch mdsw(1) 0236 case 'l' 0237 fx=linspace(log10(flmin),log10(fnyq),nfrq); % 20 Hz to Nyquist 0238 case 'm' 0239 fx=linspace(0,v_frq2mel(fnyq),nfrq); % DC to Nyquist 0240 case 'b' 0241 fx=linspace(0,v_frq2bark(fnyq),nfrq); % DC to Nyquist 0242 case 'e' 0243 fx=linspace(0,v_frq2erb(fnyq),nfrq); % DC to Nyquist 0244 otherwise % linear Hz scale 0245 fx=(0:nfrq-1)*fnyq/(nfrq-1); 0246 end 0247 else 0248 if any(mode=='h') 0249 switch mdsw(1) 0250 case 'l' 0251 fmaxu=log10(fmax); % 20 Hz to Nyquist 0252 case 'm' 0253 fmaxu=v_frq2mel(fmax); % DC to Nyquist 0254 case 'b' 0255 fmaxu=v_frq2bark(fmax); % DC to Nyquist 0256 case 'e' 0257 fmaxu=v_frq2erb(fmax); % DC to Nyquist 0258 otherwise 0259 fmaxu=fmax; % linear Hz scale 0260 end 0261 else 0262 fmaxu=fmax; % already in the correct units 0263 end 0264 if numel(fmax)<2 % only max value specified 0265 if mdsw(1)=='l' 0266 fx=linspace(log10(flmin),fmaxu,nfrq); % 20 Hz to fmax 0267 else 0268 fx=linspace(0,fmaxu,nfrq); % DC to fmax 0269 end 0270 elseif numel(fmax)<3 % min and max values specified 0271 fx=linspace(fmaxu(1),fmaxu(2),nfrq); % fmin to fmax 0272 else 0273 fmaxu(2)=fmax(2)*(fmaxu(3)-fmaxu(1))/(fmax(3)-fmax(1)); % scale the step size appropriately 0274 fx=fmaxu(1):fmaxu(2):fmaxu(3); % fmin to fmax in steps of finc 0275 nfrq=length(fx); 0276 end 0277 end 0278 switch mdsw(1) % convert the frequency range to Hz 0279 case 'l' 0280 f=10.^fx; 0281 frlab='log_{10}Hz'; 0282 frlabf='log'; 0283 frq2y=@log10; 0284 y2frq=@(x) 10.^x; 0285 case 'm' 0286 f=v_mel2frq(fx); 0287 frlab='Mel'; 0288 frlabf='Mel'; 0289 frq2y=@v_frq2mel; 0290 y2frq=@v_mel2frq; 0291 case 'b' 0292 f=v_bark2frq(fx); 0293 frlab='Bark'; 0294 frlabf='Bark'; 0295 frq2y=@v_frq2bark; 0296 y2frq=@v_bark2frq; 0297 case 'e' 0298 f=v_erb2frq(fx); 0299 frlab='Erb-rate'; 0300 frlabf='Erb'; 0301 frq2y=@v_frq2erb; 0302 y2frq=@v_erb2frq; 0303 otherwise 0304 f=fx; 0305 frlab='Hz'; 0306 frq2y=@(x) x; 0307 y2frq=@(x) x; 0308 end 0309 if ~any(mode=='H') 0310 f=fx; % give output frequencies in native units instead of Hz unless 'H' is specified 0311 end 0312 % 0313 % now calculate the spectrogram 0314 % 0315 if ns2==1 % input is a speech signal vector 0316 winlen = fix(1.81*fs(1)/bw); % window length 0317 win=0.54+0.46*cos((1-winlen:2:winlen)*pi/winlen); % Hamming window 0318 ninc=max(round(tinc(1)*fs(1)),1); % window increment in samples 0319 % we need to take account of minimum freq increment + make it exact if possible 0320 fftlen=pow2(nextpow2(4*winlen)); % enough oversampling to get good interpolation 0321 win=win/sqrt(sum(win.^2)); % ensure window squared sums to unity 0322 ix1=max(round((tinc(2)-fs(2))*fs(1)-(winlen-3)/2),1); % first sample required 0323 ix2=min(ceil((tinc(3)-fs(2))*fs(1)+(winlen+1)/2),ns1); % last sample required 0324 [sf,t]=v_enframe(s(ix1:ix2),win,ninc); 0325 t=fs(2)+(t+ix1-2)/fs(1); % time axis 0326 b=v_rfft(sf,fftlen,2); 0327 b=b.*conj(b)*2/fs(1); % Power per Hz 0328 b(:,1)=b(:,1)*0.5; % correct for no negative zero frequency to double the power 0329 b(:,end)=b(:,end)*0.5; % correct for no negative nyquist frequency to double the power 0330 fb=(0:fftlen/2)*fs(1)/fftlen; % fft bin frequencies 0331 fftfs=fs(1); 0332 else 0333 b=s; 0334 t=fs(2)+(0:ns1-1)/fs(1); % frame times 0335 fb=fs(4)+(0:ns2-1)*fs(3); 0336 fftlen=[ns2 fs(3) fs(4)]; % for v_filtbankm: ns2=# input freq bins, freq increment (Hz), first bin freq (Hz) 0337 fftfs=0; 0338 % fftlen=2*(ns2-1); % assume an even length fft 0339 % fftfs=fftlen*fs(3); 0340 end 0341 nfr=numel(t); % number of frames 0342 dblab='Power/Hz'; 0343 switch mdsw(2) 0344 case {'p','l'} 0345 b=b.*repmat(fb*log(10),nfr,1); % convert to power per decade 0346 dblab='Power/Decade'; 0347 case 'm' 0348 b=b.*repmat((700+fb)*log(1+1000/700)/1000,nfr,1); % convert to power per mel 0349 dblab='Power/Mel'; 0350 case 'b' 0351 b=b.*repmat((1960+fb).^2/52547.6,nfr,1); % convert to power per bark 0352 dblab='Power/Bark'; 0353 case 'e' 0354 b=b.*repmat(6.23*fb.^2 + 93.39*fb + 28.52,nfr,1); % convert to power per erb 0355 dblab='Power/Erb-rate'; 0356 end 0357 % 0358 % Now map onto the desired frequency scale 0359 % 0360 if any(mode=='n') 0361 fbopt=['cushn' mdsw(1)]; 0362 else 0363 fbopt=['cush' mdsw(1)]; 0364 end 0365 b=b*filtbankm(nfrq,fftlen,fftfs,fx(1),fx(end),fbopt)'; 0366 0367 if ~nargout || any(mode=='g') || any(lower(mode)=='d') 0368 if numel(db)<2 % find clipping limits 0369 plim=max(b(:))*[0.1^(0.1*db) 1]; 0370 else 0371 plim=10.^(0.1*db(1:2)); 0372 end 0373 if plim(2)<=0 0374 plim(2)=1; 0375 end 0376 if plim(1)<=0 || plim(1)==plim(2) 0377 plim(1)=0.1*plim(2); 0378 end 0379 if ~nargout || any(mode=='g') 0380 bd=10*log10(max(b,max(b(:)*1e-30))); % save an unclipped log version for plotting 0381 end 0382 if any(mode=='D') 0383 b=min(max(b,plim(1)),plim(2)); % clip the output 0384 end 0385 if any(mode=='d') 0386 b=10*log10(b); % output the dB version 0387 end 0388 end 0389 % now plot things 0390 if ~nargout || any(mode=='g') 0391 cla; % clear current axis 0392 imagesc(t,fx,bd'); 0393 axis('xy'); 0394 set(gca,'tickdir','out','clim',10*log10(plim)); 0395 if any(mode=='j') 0396 colormap('jet'); 0397 map=colormap; 0398 elseif any(mode=='J') 0399 map=tcmap; 0400 else 0401 map = repmat((0:63)'/63,1,3); 0402 end 0403 if any(mode=='i') % 'i' option = invert the colourmap 0404 map=map(64:-1:1,:); 0405 end 0406 colormap(map); 0407 if any(mode=='c') % 'c' option = show a colourbar 0408 colorbar; 0409 v_cblabel([dblab ' (dB)']); 0410 end 0411 % 0412 % Now check if annotations or a waveform are required 0413 % 0414 dotaw=[((any(mode=='t') && size(ann,2)>1) || size(ann,2)==1) size(ann,2)>1 (any(mode=='w') && ns2==1)]; 0415 ylim=get(gca,'ylim'); 0416 if any(dotaw) 0417 yrange = ylim(2)-ylim(1); 0418 zlim=ylim; 0419 toptaw=cumsum([0 dotaw.*[0.05 0.05 0.1]]*yrange)+ylim(2); 0420 zlim(2)=toptaw(4); 0421 set(gca,'ylim',zlim,'color',map(1,:)); 0422 if dotaw(3) % Plot the waveform 0423 six=min(max(floor((get(gca,'xlim')-fs(2))*fs(1))+[1 2],1),ns1); 0424 smax=max(s(six(1):six(2))); 0425 smin=min(s(six(1):six(2))); 0426 if smax==smin 0427 smax=smax+1; 0428 smin=smin-1; 0429 end 0430 srange=smax-smin; 0431 hold on 0432 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,:)) 0433 hold off 0434 end 0435 if dotaw(1) || dotaw(2) 0436 tmk=cell2mat(ann(:,1)); 0437 tmksel=tmk(:,1)<=t(end) & tmk(:,end)>=t(1); 0438 yix=1+[tmk(tmksel,1)<t(1) ones(sum(tmksel),2) tmk(tmksel,end)>t(end)]'; 0439 tmk(:,1)=max(tmk(:,1),t(1)); % clip to axis limits 0440 tmk(:,end)=min(tmk(:,end),t(end)); 0441 end 0442 if dotaw(1) && any(tmksel) % draw time markers 0443 ymk=toptaw(1:2)*[0.8 0.4;0.2 0.6]; 0444 switch size(tmk,2) 0445 case 0 0446 case 1 % isolated marks 0447 hold on 0448 plot([tmk(tmksel) tmk(tmksel)]',repmat(ymk',1,sum(tmksel)),'color',map(48,:)); 0449 hold off 0450 otherwise % draw durations 0451 0452 hold on 0453 plot(tmk(tmksel,[1 1 2 2])',ymk(yix),'color',map(48,:)); 0454 hold off 0455 end 0456 end 0457 if dotaw(2) && any(tmksel) % print annotations 0458 if any(mode=='a') 0459 horal='center'; 0460 tmk=(tmk(:,1)+tmk(:,end))*0.5; 0461 else 0462 horal='left'; 0463 tmk=tmk(:,1); 0464 end 0465 if size(ann,2)>2 0466 font='Arial'; 0467 for i=1:size(ann,1) 0468 if tmksel(i) 0469 if ~isempty(ann{i,3}) 0470 font = ann{i,3}; 0471 end 0472 text(tmk(i),toptaw(2),ann{i,2},'color',map(48,:),'fontname',font,'VerticalAlignment','baseline','HorizontalAlignment',horal); 0473 end 0474 end 0475 else 0476 for i=1:size(ann,1) 0477 if tmksel(i) 0478 text(tmk(i),toptaw(2),ann{i,2},'color',map(48,:),'VerticalAlignment','baseline','HorizontalAlignment',horal); 0479 end 0480 end 0481 end 0482 end 0483 end 0484 xlabel(['Time (' v_xticksi 's)']); 0485 if any(mode=='f') && ~strcmp(frlab,'Hz') 0486 ylabel([frlabf '-scaled frequency (Hz)']); 0487 ytickhz(frq2y,y2frq); 0488 else 0489 ylabel(['Frequency (' v_yticksi frlab ')']); 0490 end 0491 ytick=get(gca,'YTick'); 0492 ytickl=get(gca,'YTickLabel'); 0493 msk=ytick<=ylim(2); 0494 if any(~msk) 0495 set(gca,'YTick',ytick(msk),'YTickLabel',ytickl(msk)); 0496 end 0497 end 0498 0499 function ytickhz(frq2y,y2frq) 0500 % label non linear y frequency axis 0501 % 0502 % Bugs/Suggestions: 0503 % * Add a penalty for large numbers (e.g. 94 is less "round" than 11) 0504 % * possibly add subticks at 1:2:5 if boundaries are 1 and 10 0505 % * could treat subtick allocation specially if bounding lables are both powers of 10 0506 % and work in log spacing rather than spacing directly 0507 0508 % algorithm constants 0509 0510 seps=[0.4 1 3 6]; % spacings: (a) min subtick, (b) min tick, (c) min good tick, (d) max good tick 0511 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) 0512 nbest=10; % number of possibilities to track 0513 0514 prefix={'y','z','a','f','p','n','u','m','','k','M','G','T','P','E','Z','Y'}; 0515 0516 ah=gca; 0517 getgca=get(ah); % Get original axis properties 0518 set(ah,'Units','points','FontUnits','points'); 0519 getgcac=get(ah); % Get axis properties in points units 0520 set(ah,'Units',getgca.Units,'FontUnits',getgca.FontUnits); % return to original values 0521 ylim=getgca.YLim; 0522 yrange=ylim*[-1;1]; 0523 chsz= yrange*getgcac.FontSize/getgcac.Position(4); % char height in Y-units 0524 % divide the y-axis up into bins containing at most one label each 0525 maxl=ceil(2*yrange/chsz); % max number of labels 0526 0527 % candidate array [cand(:,[1 2])/1000 cand(:,5) cand(:,6)/1000 cand(:,[7 8])] 0528 % 1,2=y limits, 3,4=log limits, 5=Hz, 6=cost, 7=mantissa, 8=exponent, 9=sig digits, 10=y-position 0529 cand=zeros(maxl+2,10); 0530 yinc=(yrange+chsz*0.0002)/maxl; % bin spacing (allowing for a tiny bit to ensure the ends are included) 0531 cand(2:end-1,2)=ylim(1)+yinc*(1:maxl)'-chsz*0.0001; 0532 cand(3:end-1,1)=cand(2:end-2,2); 0533 cand(2,1)=cand(2,2)-yinc; 0534 cand(2:end-1,1:2)=y2frq(max(cand(2:end-1,1:2),0)); 0535 0536 % find the "roundest" number in each interval 0537 % first deal with intervals containing zero 0538 cand([1 maxl+2],6)=-1; 0539 cand(2,9)=(cand(2,1)<=0); % mask out interval contaiing zero 0540 cand(2,6)=-cand(2,9); 0541 msk=cand(:,6)==0; % find rows without a cost yet 0542 cand(msk,3:4)=log10(cand(msk,1:2)); 0543 % find powers of 1000 0544 loglim=ceil(cand(:,3:4)/3); 0545 msk=loglim(:,2)>loglim(:,1); 0546 if any(msk) 0547 xp=loglim(msk,1); 0548 wuns=ones(length(xp),1); 0549 cand(msk,5:9)=[1000.^xp wuns-ww(3) wuns 3*xp wuns]; 0550 end 0551 % find powers of 10 0552 loglim=ceil(cand(:,3:4)); 0553 msk=~msk & (loglim(:,2)>loglim(:,1)); 0554 if any(msk) 0555 xp=loglim(msk,1); 0556 wuns=ones(length(xp),1); 0557 cand(msk,5:9)=[10.^xp wuns-ww(2) wuns xp wuns]; 0558 end 0559 % find value with fewest digits 0560 msk=cand(:,6)==0; % find rows without a cost yet 0561 maxsig=1-floor(log10(10^min(cand(msk,3:4)*[-1;1])-1)); % maximum number of significant figures to consider 0562 pten=10.^(0:maxsig-1); % row vector of powers of ten 0563 noten=10.^(-floor(cand(msk,3))); % exponent of floating point representation of lower bound 0564 sigdig=sum((ceil(cand(msk,2).*noten*pten)-ceil(cand(msk,1).*noten*pten))==0,2); % number of digits common to the interval bounds 0565 lowman=ceil(cand(msk,1).*noten.*10.^sigdig); 0566 midman=10*floor(lowman/10)+5; 0567 highman=ceil(cand(msk,2).*noten.*10.^sigdig); 0568 mskman=midman>=lowman & midman<highman; % check if we can include a manitssa ending in 5 0569 lowman(mskman)=midman(mskman); 0570 cand(msk,6:9)=[sigdig+1 lowman floor(cand(msk,3))-sigdig sigdig+1]; 0571 cand(msk,5)=cand(msk,7).*10.^cand(msk,8); 0572 cand(msk,6)=cand(msk,6)-(mod(cand(msk,7),10)==5)*ww(1); 0573 cand(2:end-1,10)=frq2y(cand(2:end-1,5)); 0574 cand([1 maxl+2],10)=ylim + seps(4)*chsz*[-1 1]; % put imaginary labels at the optimum spacing beyond the axes 0575 % [cand(:,[1 2 5])/1000 cand(:,[6 7 8 9])] 0576 0577 % Now do n-best DP to find the best sequence 0578 0579 ratint=[8/5 25/10 0 0 4/3]; 0580 costs=Inf(nbest,maxl+2); % cumulative path costs 0581 costs(1,1)=0; % starting node only has one option 0582 prev=ones(nbest,maxl+2); % previous label in path 0583 labcnt=zeros(nbest,maxl+2); % number of labels in path 0584 for i=2:maxl+2 0585 ntry=nbest*(i-1); % number of previous options 0586 prevc=reshape(repmat(1:i-1,nbest,1),ntry,1); % previous candidate 0587 prevprev=1+floor((prev(1:ntry)'-1)/nbest); % previous previous candidate 0588 msk=prevprev>1+(maxl+2)*(i==maxl+2); % mask for label triplets 0589 labcnti=labcnt(1:ntry)+1; 0590 disti=(cand(i,10)-cand(prevc,10))/chsz; % distance to previous label in characters 0591 costa=max(seps(3)-disti,0)*ww(6)+max(disti-seps(4),0)*ww(7); 0592 incri=(cand(i,5)-cand(prevc,5)); % label increment 0593 incrj=(cand(i,5)-cand(prevprev,5)); % double label increment 0594 if any(msk) 0595 costa(msk)=costa(msk)- ww(4)*(abs(incrj(msk)-2*incri(msk))<0.01*incri(msk)); 0596 if cand(i,7)==1 || cand(i,7)==2 || cand(i,7)==5 % look for labels 1:2:5 0597 costa(msk)=costa(msk)- ww(5)*(abs(incrj(msk)-ratint(cand(i,7))*incri(msk))<0.01*incri(msk)); 0598 end 0599 end 0600 costa(disti<seps(2))=Inf; 0601 costi=(costs(1:ntry).*max(labcnt(1:ntry),1)+costa'+cand(i,6))./labcnti; 0602 [sc,isc]=sort(costi); 0603 isc=isc(1:nbest); 0604 costs(:,i)=sc(1:nbest)'; 0605 prev(:,i)=isc'; 0606 labcnt(:,i)=labcnti(isc)'; 0607 end 0608 0609 % now traceback the best sequence 0610 0611 % fprintf('Traceback\n\n'); 0612 ichoose=0; 0613 labchoose=[]; 0614 for i=1:nbest 0615 if labcnt(i,maxl+2)>1 && costs(i,maxl+2)<Inf 0616 lablist=zeros(labcnt(i,maxl+2)-1,1); 0617 k=prev(i,maxl+2); 0618 for j=labcnt(i,maxl+2)-1:-1:1 0619 lablist(j)=1+floor((k-1)/nbest); 0620 k=prev(k); 0621 end 0622 % fprintf('Cost=%8.2f :',costs(i,maxl+2)); 0623 % fprintf(' %g',cand(lablist,5)) 0624 % fprintf('\n'); 0625 if ~ichoose || labcnt(ichoose,maxl+2)==1 0626 ichoose=i; 0627 labchoose=lablist; 0628 end 0629 end 0630 end 0631 0632 % now create the labels 0633 0634 ntick=length(labchoose); 0635 % sort out the subticks 0636 subpos=[]; 0637 if ntick>=2 0638 for i=1:ntick-1 0639 clj=cand(labchoose(i:i+1),:); 0640 sprec=min(clj(1,8)+100*(clj(1,7)==0),clj(2,8)); % subtick precision 0641 spos=(clj(1,7)*10^(clj(1,8)-sprec):clj(2,7)*10^(clj(2,8)-sprec))*10^sprec; 0642 nsub=length(spos); 0643 if nsub==2 0644 spos=spos*[1 0.5 0;0 0.5 1]; 0645 nsub=3; 0646 end 0647 if nsub>=3 0648 yspos=frq2y(spos); 0649 for kk=1:3 % try various subdivisions: every 1, 2 or 5 0650 k=kk+2*(kk==3); % 1, 2 and 5 0651 if 2*k<=nsub-1 && ~mod(nsub-1,k) % must divide exactly into nsub 0652 if all((yspos(1+k:k:nsub)-yspos(1:k:nsub-k))>=(seps(1)*chsz)) % check they all fit in 0653 subpos=[subpos yspos(1+k:k:nsub-k)]; 0654 if i==1 0655 spos=(ceil(cand(2,1)/10^sprec):clj(1,7)*10^(clj(1,8)-sprec))*10^sprec; 0656 nsub=length(spos); 0657 yspos=frq2y(spos); 0658 if nsub>=k+1 && all((yspos(nsub:-k:1+k)-yspos(nsub-k:-k:1))>=(seps(1)*chsz)) 0659 subpos=[subpos yspos(nsub-k:-k:1)]; 0660 end 0661 elseif i==ntick-1 0662 spos=(clj(2,7)*10^(clj(2,8)-sprec):floor(cand(end-1,2)/10^sprec))*10^sprec; 0663 nsub=length(spos); 0664 yspos=frq2y(spos); 0665 if nsub>=k+1 && all((yspos(1+k:k:nsub)-yspos(1:k:nsub-k))>=(seps(1)*chsz)) 0666 subpos=[subpos yspos(1+k:k:nsub)]; 0667 end 0668 end 0669 break; 0670 end 0671 end 0672 end 0673 end 0674 end 0675 end 0676 nsub=length(subpos); 0677 tickpos=[cand(labchoose,10); subpos']; 0678 ticklab=cell(ntick+nsub,1); 0679 sipref=min(max(floor((sum(cand(labchoose,8:9),2)-1)/3),-8),8); 0680 nzadd=cand(labchoose,8)-3*sipref; % trailing zeros to add 0681 digzer=cand(labchoose,7).*10.^max(nzadd,0); % label digits including trailing zeros 0682 ndleft=cand(labchoose,9)+nzadd; % digits to the left of the decimal point 0683 for i=1:ntick 0684 tickint=num2str(digzer(i)); 0685 if nzadd(i)<0 0686 tickint=[tickint(1:ndleft(i)) '.' tickint(1+ndleft(i):end)]; 0687 end 0688 ticklab{i} = sprintf('%s%s',tickint,prefix{sipref(i)+9}); 0689 end 0690 for i=ntick+1:ntick+nsub 0691 ticklab{i}=''; 0692 end 0693 [tickpos,ix]=sort(tickpos); 0694 ticklab=ticklab(ix); 0695 0696 set(ah,'YTick',tickpos','YTickLabel',ticklab); 0697