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