STDSPECTRUM Generate standard acoustic/speech spectra in s- or z-domain [B,A]=(S,M,F,N,ZI,BS,AS)
Inputs: s Spectrum type (either text or number - see below) or 0 to use bs/as
m mode: char 1 specifies output type,
f - frequency response (complex)
m - magnitude response
p - power spectrum
d - decibel power spectrum
t - time waveform
s - s-domain filter: b(s)/a(s) [default]
z - z-domain filter: b(z)/a(z)
i - sampled impulse response
plotting options
M - plot magnitude spectrum
E - plot magnitude spectrum error
Q - plot phase spectrum
T - also plot target spectra
A - plot zeros/poles
U - plot using a uniform frequency axis
L - plot using a log frequency axis
W - waveform
S - spectrogram
f sample frequency (modes z,i,t) or set of frequencies in Hz (modes f,m,p,d)
n number of output samples (mode i,t)
zi initial state of filter from a previous call (mode t)
bs numerator s-domain polynomial (or cell array containing polynomial factors)
as denominator s-domain polynomial (or cell array containing polynomial factors)
Outputs: The outputs depend on the mode selected:
mode = 'f', 'm', 'p' or 'd'
b = ouptut spectrum
mode = 's' or 'z'
b,a = numerator and denonminator of the output spectrum
mode = 't'
b = output waveform
a = final state of the filter - use as the zi input of a future call
Spectrum types (specify either as a number or case-insensitive text abbreviation):
0 external : BS and AS arguments specify an s-domain filter
1 White : white noise
2 A-Weight : the formula for this is given in [3] and is based on
the equal-loudness curves of [9]
3 B-Weight : this needs to be confirmed with ANSI S1.4-1981 standard or IEC 60651
4 C-Weight : the formula for this is given in [3]
5 LTASS-P50 : the long-term average speech spectrum that is defined by a
formula on page 3 of [4] which, strangely, does not precisely
match the graph shown on the same page.
6 LTASS-1994 : the long-term average speech spectrum that is taken from a graph in [2]
7 SII-intinv : The inverse spectrum of the ear's internal masking noise; this is taken
from table 1 of [1]. It is inverted so that it is a bandpass rather than
bandstop characteristic.
8 BS-468 : The weighting proposed for audio frequency noise measurement in [5] and [6].
9 USASI : Noise simulating long-term programme material spectrum from [7],[8].
The level is such that the power is 0dB over an infinite bandwidth
10 POTS : the D spectrum from [11].
References:
[1] Methods for the calculation of the speech intelligibility index.
ANSI Standard S3.5-1997 (R2007), American National Standards Institute, 1997.
[2] D. Byrne, H. Dillon, K. Tran, S. Arlinger, K. Wilbraham, R. Cox, B. Hayerman,
R. Hetu, J. Kei, C. Lui, J. Kiessling, M. N. Kotby, N. H. A. Nasser,
W. A. H. E. Kholy, Y. Nakanishi, H. Oyer, R. Powell, D. Stephens, R. Meredith,
T. Sirimanna, G. Tavartkiladze, G. I. Frolenkov, S. Westerman, and C. Ludvigsen.
An international comparison of long-term average speech spectra.
JASA, 96 (4): 2108–2120, Oct. 1994.
[3] CENELEC. Electroacoustics - sound level meters. Technical Report EN EN 61672-1:2003, 2003.
(also ANSI S1.42-2001)
[4] ITU-T. Artificial voices. Standard P.50, Sept. 1999.
[5] ITU-T. Measurement of weighted noise in sound-programme circuits.
Recommendation J.16, 1988.
[6] ITU-R. Measurement of audio-requency noise voltage level in sound broadcasting.
Recommendation BS.468., 1986.
[7] NRSC AM Reemphasis, Deemphasize, and Broadcast Audio Transmission Bandwidth Specifications,
EIA-549 Standard, Electronics Industries Association , July 1988.
[8] NRSC AM Reemphasis, Deemphasize, and Broadcast Audio Transmission Bandwidth Specifications,
NRSC-1-A Standard, Sept 2007, Online: http://www.nrscstandards.org/SG/NRSC-1-A.pdf
[9] H. Fletcher and W. A. Munson. Loudness, its definition, measurement and calculation.
J. Acoust Soc Amer, 5: 82–108, Oct. 1933.
[10] American National Standard Specification for Sound Level Meters.
ANSI S1.4-1983 (R2006)/ANSI S1.4a-1985 (R2006), American National Standards Institute
[11] IEEE standard equipment requirements and measurement techniques for analog transmission
parameters for telecommunications. Standard IEEE Std 743-1995, Dec. 1995.0001 function [b,a]=stdspectrum(s,m,f,n,zi,bs,as) 0002 %STDSPECTRUM Generate standard acoustic/speech spectra in s- or z-domain [B,A]=(S,M,F,N,ZI,BS,AS) 0003 % 0004 %Inputs: s Spectrum type (either text or number - see below) or 0 to use bs/as 0005 % m mode: char 1 specifies output type, 0006 % f - frequency response (complex) 0007 % m - magnitude response 0008 % p - power spectrum 0009 % d - decibel power spectrum 0010 % t - time waveform 0011 % s - s-domain filter: b(s)/a(s) [default] 0012 % z - z-domain filter: b(z)/a(z) 0013 % i - sampled impulse response 0014 % plotting options 0015 % M - plot magnitude spectrum 0016 % E - plot magnitude spectrum error 0017 % Q - plot phase spectrum 0018 % T - also plot target spectra 0019 % A - plot zeros/poles 0020 % U - plot using a uniform frequency axis 0021 % L - plot using a log frequency axis 0022 % W - waveform 0023 % S - spectrogram 0024 % f sample frequency (modes z,i,t) or set of frequencies in Hz (modes f,m,p,d) 0025 % n number of output samples (mode i,t) 0026 % zi initial state of filter from a previous call (mode t) 0027 % bs numerator s-domain polynomial (or cell array containing polynomial factors) 0028 % as denominator s-domain polynomial (or cell array containing polynomial factors) 0029 % 0030 % Outputs: The outputs depend on the mode selected: 0031 % mode = 'f', 'm', 'p' or 'd' 0032 % b = ouptut spectrum 0033 % mode = 's' or 'z' 0034 % b,a = numerator and denonminator of the output spectrum 0035 % mode = 't' 0036 % b = output waveform 0037 % a = final state of the filter - use as the zi input of a future call 0038 % 0039 % Spectrum types (specify either as a number or case-insensitive text abbreviation): 0040 % 0 external : BS and AS arguments specify an s-domain filter 0041 % 1 White : white noise 0042 % 2 A-Weight : the formula for this is given in [3] and is based on 0043 % the equal-loudness curves of [9] 0044 % 3 B-Weight : this needs to be confirmed with ANSI S1.4-1981 standard or IEC 60651 0045 % 4 C-Weight : the formula for this is given in [3] 0046 % 5 LTASS-P50 : the long-term average speech spectrum that is defined by a 0047 % formula on page 3 of [4] which, strangely, does not precisely 0048 % match the graph shown on the same page. 0049 % 6 LTASS-1994 : the long-term average speech spectrum that is taken from a graph in [2] 0050 % 7 SII-intinv : The inverse spectrum of the ear's internal masking noise; this is taken 0051 % from table 1 of [1]. It is inverted so that it is a bandpass rather than 0052 % bandstop characteristic. 0053 % 8 BS-468 : The weighting proposed for audio frequency noise measurement in [5] and [6]. 0054 % 9 USASI : Noise simulating long-term programme material spectrum from [7],[8]. 0055 % The level is such that the power is 0dB over an infinite bandwidth 0056 % 10 POTS : the D spectrum from [11]. 0057 % 0058 % References: 0059 % [1] Methods for the calculation of the speech intelligibility index. 0060 % ANSI Standard S3.5-1997 (R2007), American National Standards Institute, 1997. 0061 % [2] D. Byrne, H. Dillon, K. Tran, S. Arlinger, K. Wilbraham, R. Cox, B. Hayerman, 0062 % R. Hetu, J. Kei, C. Lui, J. Kiessling, M. N. Kotby, N. H. A. Nasser, 0063 % W. A. H. E. Kholy, Y. Nakanishi, H. Oyer, R. Powell, D. Stephens, R. Meredith, 0064 % T. Sirimanna, G. Tavartkiladze, G. I. Frolenkov, S. Westerman, and C. Ludvigsen. 0065 % An international comparison of long-term average speech spectra. 0066 % JASA, 96 (4): 2108–2120, Oct. 1994. 0067 % [3] CENELEC. Electroacoustics - sound level meters. Technical Report EN EN 61672-1:2003, 2003. 0068 % (also ANSI S1.42-2001) 0069 % [4] ITU-T. Artificial voices. Standard P.50, Sept. 1999. 0070 % [5] ITU-T. Measurement of weighted noise in sound-programme circuits. 0071 % Recommendation J.16, 1988. 0072 % [6] ITU-R. Measurement of audio-requency noise voltage level in sound broadcasting. 0073 % Recommendation BS.468., 1986. 0074 % [7] NRSC AM Reemphasis, Deemphasize, and Broadcast Audio Transmission Bandwidth Specifications, 0075 % EIA-549 Standard, Electronics Industries Association , July 1988. 0076 % [8] NRSC AM Reemphasis, Deemphasize, and Broadcast Audio Transmission Bandwidth Specifications, 0077 % NRSC-1-A Standard, Sept 2007, Online: http://www.nrscstandards.org/SG/NRSC-1-A.pdf 0078 % [9] H. Fletcher and W. A. Munson. Loudness, its definition, measurement and calculation. 0079 % J. Acoust Soc Amer, 5: 82–108, Oct. 1933. 0080 % [10] American National Standard Specification for Sound Level Meters. 0081 % ANSI S1.4-1983 (R2006)/ANSI S1.4a-1985 (R2006), American National Standards Institute 0082 % [11] IEEE standard equipment requirements and measurement techniques for analog transmission 0083 % parameters for telecommunications. Standard IEEE Std 743-1995, Dec. 1995. 0084 0085 % Other candidates: (a) Z-weighting, (b) ISO226, (c) USASI, (d) P.48 spectra 0086 % 0087 % Other standards: 0088 % IEEE743 has several weighting filters defined 0089 % ITU-T 0.41 Psophometer for use on telephone-type circuits 0090 % Bell System Technical Reference 41009 (C-message) 0091 % ISO 8041:2005 (E): Human Response to Vibration – Measuring 0092 % Instrumentation 0093 % IEC 1260:1995, class 1 (also IEC 61260/ANSI S1.11-2004) Octave band and fractional octave band filters 0094 % IEC 651: Specification for Sound Level Meters 0095 % IRS P.48 0096 % mIRS P.830 0097 0098 % Copyright (C) Mike Brookes 2008 0099 % Version: $Id: stdspectrum.m 1243 2011-12-28 11:55:53Z dmb $ 0100 % 0101 % VOICEBOX is a MATLAB toolbox for speech processing. 0102 % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 0103 % 0104 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0105 % This program is free software; you can redistribute it and/or modify 0106 % it under the terms of the GNU General Public License as published by 0107 % the Free Software Foundation; either version 2 of the License, or 0108 % (at your option) any later version. 0109 % 0110 % This program is distributed in the hope that it will be useful, 0111 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0112 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0113 % GNU General Public License for more details. 0114 % 0115 % You can obtain a copy of the GNU General Public License from 0116 % http://www.gnu.org/copyleft/gpl.html or by writing to 0117 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0118 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0119 0120 % Bugs/Suggestions: 0121 % * Could generate impulse response by appending a LP filter (elliptic) to 0122 % the s-domain transfer function and sampling the impulse response 0123 % * better calculation of impulse response length based on its total power 0124 % * ensure that the number of z-domain zeros at z=1 is correct. 0125 0126 persistent spty spz az bz fz ixsz maxj 0127 % spty contains the name of the spectrum 0128 % spz contains a list of the poles and zeros 0129 % spz(1) = gain, spz(2) = number of zeros (excluding implicit conjugates), spz(3...) zeros followed by poles 0130 % only one of a complex conjugate pair is given 0131 if isempty(spz) 0132 spty={'White';'A-Weight';'B-Weight';'C-Weight';'LTASS-P50';'LTASS-1994';'SII-IntInv';'BS-468';'USASI';'POTS'}; 0133 spz={[1 0]; 0134 [7390100803.6603 4 0 0 0 0 -129.42731565506 -129.42731565506 -676.40154023295 -4636.125126885 -76618.526016858 -76618.526016858]; 0135 [5986190155.0156 3 0 0 0 -129.42731565506 -129.42731565506 -995.88487118796 -76618.526016858 -76618.526016858]; 0136 [5912384617.784 2 0 0 -129.42731565506 -129.42731565506 -76618.526016858 -76618.526016858]; 0137 [1.1294790345421e+015 3 0 0 -34437.856184098 -721.94747118664+754.97798119504i -1721.704402273 -5234.2950286868 -10953.570773216+42789.342252749i]; 0138 [19720493.192959 5 0 0 0 -22550.637578954 -11319.635610404+70239.177107659i -253.31327846696+672.10855509952i -1299.1885437118+2301.2064056419i -10646.952627978+68290.702816027i -147307.51763333]; 0139 [6.1311266394354e+018 2 0 0 -381.08293630892 -5920.974797779 -4701.76218192+24369.279310049i 10597.854874768+39258.915154617i]; 0140 [2.1034520039796e+024 1 0 -25903.701047817 -23615.535213635+36379.908937329i -62675.170058468 -18743.746690721+62460.156452506i]; 0141 [72.648989380657 1 0 -2*pi*[100 320]]; 0142 [7.8820088171767e+016 4 0 0 0 0 -452.681+1924.28i -2334+1702.73i -11264.2+8213.32i -4665.8+19828.7i]; 0143 }; 0144 fz=-1; 0145 ixsz=-1; 0146 end 0147 if nargin<2 || ~numel(m) 0148 if nargout 0149 m='s '; % default is Laplace transform with default unit 0150 else 0151 m='d '; 0152 end 0153 end 0154 m1=m(1); % output format 0155 if nargin<3 0156 f=8192; % default frequency 0157 end 0158 0159 % determine the spectrum type 0160 0161 if ~numel(s) || s(1)==0 0162 ixs=0; 0163 sb=1; 0164 sz=[]; % list of s-domain zeros 0165 if iscell(bs) 0166 for i=1:numel(bs) 0167 sb=conv(sb,bs{i}); 0168 sz=[sz roots(bs{i}).']; 0169 end 0170 elseif numel(bs)>0 0171 sb=bs; 0172 sz=roots(bs).'; 0173 end 0174 sa=1; 0175 sp=[]; % list of s-domain poles 0176 if iscell(as) 0177 for i=1:numel(as) 0178 sa=conv(sa,as{i}); 0179 sp=[sp roots(as{i}).']; 0180 end 0181 elseif numel(as)>0 0182 sa=as; 0183 sp=roots(as).'; 0184 end 0185 else 0186 if ischar(s) 0187 ixs=find(strcmpi(s,spty)); 0188 if isempty(ixs) 0189 error('undefined spectrum type: %s',s); 0190 end 0191 else 0192 ixs=s; 0193 end 0194 if ixs>size(spty,1) 0195 error('undefined spectrum type: %d',ixs); 0196 end 0197 0198 % get s-domain function 0199 % sb/sa is transfer function 0200 % sz and sp are lists of zeros and poles of lengths nsp and nsz 0201 0202 spzi=spz{ixs}; 0203 nsz=spzi(2); 0204 sz=spzi(3:3+nsz-1); 0205 sz=[sz conj(sz(imag(sz)~=0))]; 0206 sp=spzi(3+nsz:end); 0207 sp=[sp conj(sp(imag(sp)~=0))]; 0208 sb=spzi(1)*poly(sz); 0209 sa=poly(sp); 0210 end 0211 if (nargin<3 || ~numel(f)) && any(m1=='fmpd') % calcualte the frequency range 0212 apz=abs([sp sz]); 0213 apz(apz==0)=[]; % ignore zero frequency poles/zeros 0214 if ~numel(apz) 0215 apz=[100 5000]; 0216 elseif length(apz)==1 0217 apz=[apz/10 apz*10]; 0218 end 0219 f=logspace(log10(min(apz)*0.5/pi)-0.5,log10(max(apz)*0.5/pi)+0.5); 0220 end 0221 if any(m1=='fmpd') 0222 h=polyval(sb,2i*pi*f)./polyval(sa,2i*pi*f); 0223 end 0224 fs=f; % save sampling frequency 0225 if any(m1=='izt') && (f~=fz || ixs~=ixsz || ixs==0) 0226 % we use an iterative method to find the best digital filter 0227 % we initialize the phases with either bilinear or impulse invariance 0228 % only using the impulse invariance if it is good (max error < 10dB) 0229 % we then iterate invfreqz using the s-domain magnitudes and 0230 % the phases of the best fit so far. 0231 % we use log-spaced frequencies at low frequencies and linear at high 0232 % we then search for various numbers of poles and zeros 0233 fz=f; % save sampling frequency for future call 0234 ixsz=ixs; 0235 if ixs==1 0236 bz=1; 0237 az=1; 0238 maxj=0; 0239 else 0240 warning off all % avoid lots of ill-conditioning error messages 0241 nflin=100; % number of frequency samples in linear region (high freq) 0242 alp=1.15; % freq ratio increment in low freq region 0243 f0=25*2*pi/f; % minimum interesting frequency (25 Hz in radians) 0244 fx=pi/nflin/(alp-1); % boundary between log and linear portions 0245 if fx<=f0 || f0>=pi 0246 fif=linspace(0,pi,nflin); 0247 elseif fx>pi 0248 fif=[0 logspace(log10(f0),log10(pi),ceil(log10(pi/f0)/log10(alp)))]; 0249 else 0250 nlin=ceil((pi-fx)*nflin/pi); 0251 fif=[0 logspace(log10(f0),log10(fx),ceil(log10(fx/f0)/log10(alp))) linspace(fx+(pi-fx)/nlin,pi,nlin-1)]; 0252 end 0253 h0=abs(polyval(sb,1i*fif*f)./polyval(sa,1i*fif*f)); % target magnitude spectrum 0254 hix=h0~=0; % don't calculate dB errors if zero (e.g. at DC) 0255 % initialize with impulse invariance 0256 [bz,az]=impinvar(sb,sa,f); 0257 hj=freqz(bz,az,fif); 0258 maxj=max(abs(db(abs(hj(hix)))-db(abs(h0(hix))))); 0259 % or else with bilinear 0260 [ifb,ifa]=bilinear(sb,sa,f); 0261 hn=freqz(ifb,ifa,fif); 0262 maxi=max(abs(db(abs(hn(hix)))-db(abs(h0(hix))))); 0263 if maxi<maxj || maxj>10 % accept bilinear if it is better or if imp inv is bad 0264 maxj=maxi; 0265 bz=ifb; 0266 az=ifa; 0267 hj=hn; 0268 end 0269 for mm=1:length(sa) % maximum number of poles 0270 for nn=1:mm % number of zeros is always less than number of poles 0271 hn=hj; 0272 j=0; 0273 for i=1:30 % iterate up to 30 times (usually less) 0274 h=h0.*exp(1i*angle(hn)); 0275 [ifb,ifa]=invfreqz(h,fif,nn,mm,[],10); 0276 hn=freqz(ifb,ifa,fif); 0277 maxi=max(abs(db(abs(hn(hix)))-db(abs(h0(hix))))); 0278 if maxi<maxj 0279 maxj=maxi; 0280 bz=ifb; 0281 az=ifa; 0282 hj=hn; 0283 j=i; 0284 end 0285 if i>j+5 % quit if no improvement in last five iterations 0286 break 0287 end 0288 end 0289 end 0290 end 0291 warning on all 0292 end 0293 end 0294 switch m1 0295 case 'z' 0296 b=bz; 0297 a=az; 0298 case 't' 0299 if nargin<5 || ~numel(zi) 0300 [b,a]=randfilt(bz,az,n); 0301 else 0302 [b,a]=randfilt(bz,az,n,zi); 0303 end 0304 case 'i' 0305 if nargin<5 || ~numel(zi) 0306 if nargin<4 || ~numel(n) % determine n to include 1 - 1e-8 of the energy 0307 n=ceil(-fs*log(1e4)/max(real(sp))); 0308 end 0309 [b,a]=filter(bz,az,[1; zeros(n-1,1)]); 0310 else 0311 [b,a]=filter(bz,az,zeros(n,1),zi); 0312 end 0313 case 'm' 0314 b = abs(h); 0315 case 'f' 0316 b = h; 0317 case 'd' 0318 b = db(abs(h)); 0319 case 'p' 0320 b=h.*conj(h); 0321 case 's' 0322 b=sb; 0323 a=sa; 0324 otherwise 0325 error('Output format %s not yet implemented',m1); 0326 end 0327 0328 % plot data 0329 if ~nargout || ~strcmp(m,lower(m)) 0330 if strcmp(m,lower(m)) % if no upper case letters 0331 m='ML'; 0332 end 0333 if ~any(m=='Q') && ~any(m=='A') && ~any(m=='E') && (m1~='t' || ~any(m=='W') && ~any(m=='S')) 0334 m(end+1)='M'; % default plot type 0335 end 0336 nfig=0; 0337 paz=any(m1=='itz'); % plot discrete time result 0338 pas=~paz || any(m=='T'); % plot continuous time result 0339 if any(m=='M') || any(m=='Q') || any(m=='E') % draw a frequency response plot 0340 clf; 0341 nfig=1; 0342 pam=any(m=='M'); % magnitude response 0343 pae=any(m=='E') && paz; % magniture response error 0344 paq=any(m=='Q'); % phase response 0345 pat=any(m=='T'); % include target spectrum 0346 pal=any(m=='L') || (~paz && ~any(m=='U')); % log frequency axis 0347 if any(m1=='itz') 0348 fs=f; % save the sample frequency 0349 apz=abs([sp sz]); % and determine the frequency range 0350 apz(apz==0)=[]; % ignore zero frequency poles/zeros 0351 if ~numel(apz) 0352 apz=[100 5000]; 0353 elseif length(apz)==1 0354 apz=[apz/10 apz*10]; 0355 end 0356 if pal 0357 f=logspace(log10(min([fs/1000 apz*0.05/pi])),log10(fs/2),200); 0358 else 0359 f=linspace(min([fs/1000 apz*0.05/pi]),fs/2,200); 0360 end 0361 end 0362 hs=freqs(sb,sa,2*pi*f); 0363 if paz 0364 hz=freqz(bz,az,f,fs); 0365 end 0366 0367 axh=[]; 0368 nax=pam+pae+paq; % number of axis sets 0369 0370 if pam 0371 if nax>1 0372 subplot(nax,1,1); 0373 end 0374 if paz 0375 if pas 0376 plot(f,db(abs(hs)),':r',f,db(abs(hz)),'-b') 0377 ymax=max(db(abs([hs hz])))+1; 0378 ymin=min(db(abs([hs hz])))-1; 0379 else 0380 plot(f,db(abs(hz)),'-b') 0381 ymax=max(db(abs(hz)))+1; 0382 ymin=min(db(abs(hz)))-1; 0383 end 0384 else 0385 plot(f,db(abs(hs)),'-b') 0386 ymax=max(db(abs(hs)))+1; 0387 ymin=min(db(abs(hs)))-1; 0388 end 0389 if pal 0390 set(gca,'Xscale','log'); 0391 end 0392 set(gca,'XLim',[min(f) max(f)],'YLim',[max(ymin,ymax-60) ymax]); 0393 xlabel(['Frequency (' xticksi 'Hz)']); 0394 ylabel('Gain (dB)'); 0395 if ixs>0 0396 title(spty{ixs}); 0397 end 0398 axh(end+1)=gca; 0399 end 0400 if pae 0401 if nax>1 0402 subplot(nax,1,1+pam); 0403 end 0404 plot(f,db(abs(hz))-db(abs(hs)),'-b') 0405 if pal 0406 set(gca,'Xscale','log'); 0407 end 0408 set(gca,'XLim',[min(f) max(f)]); 0409 xlabel(['Frequency (' xticksi 'Hz)']); 0410 ylabel('Gain Error (dB)'); 0411 if ixs>0 && ~pam 0412 title(spty{ixs}); 0413 end 0414 axh(end+1)=gca; 0415 end 0416 if paq 0417 if nax>1 0418 subplot(nax,1,nax); 0419 end 0420 if paz 0421 if pas 0422 plot(f,angle(hs),':r',f,angle(hz),'-b') 0423 else 0424 plot(f,angle(hz),'-b') 0425 end 0426 else 0427 plot(f,angle(hs),'-b') 0428 end 0429 if pal 0430 set(gca,'Xscale','log'); 0431 end 0432 set(gca,'XLim',[min(f) max(f)]); 0433 xlabel(['Frequency (' xticksi 'Hz)']); 0434 ylabel('Phase (rad)'); 0435 if ixs>0 && nax==1 0436 title(spty{ixs}); 0437 end 0438 axh(end+1)=gca; 0439 end 0440 if nax>1 0441 linkaxes(axh,'x'); 0442 end 0443 end 0444 if any(m=='A') % plot complex plane 0445 if nfig 0446 figure(); 0447 end 0448 clf; 0449 nfig=1; 0450 if pas 0451 if paz 0452 subplot(121); 0453 end 0454 plot(real(sp),imag(sp),'xb',real(sz),imag(sz),'ob'); 0455 axis equal; 0456 xlim=get(gca,'xlim'); 0457 xlim(1)=min(xlim(1),-1000); 0458 xlim(2)=max(xlim(2),1000); 0459 ylim=get(gca,'ylim'); 0460 ylim(1)=min(ylim(1),-1000); 0461 ylim(2)=max(ylim(2),1000); 0462 axis([xlim ylim]); 0463 hold on 0464 plot(xlim,[0 0],':r',[0 0],ylim,':r'); 0465 hold off 0466 title('s-plane'); 0467 end 0468 if paz 0469 if pas 0470 subplot(122); 0471 end 0472 axl=max(abs([1.1; az(:);bz(:)])); 0473 t=linspace(0,2*pi); 0474 rtzb=roots(bz); 0475 rtza=roots(az); 0476 plot(cos(t),sin(t),':r',[-1 0; 1 0],[0 -1; 0 1],':r',real(rtza),imag(rtza),'xb',real(rtzb),imag(rtzb),'ob'); 0477 axis equal; 0478 axis([-1 1 -1 1]*axl); 0479 title('z-plane'); 0480 end 0481 end 0482 if any(m=='W') && any(m1=='it') % plot waveform 0483 if nfig 0484 figure(); 0485 end 0486 clf; 0487 nfig=1; 0488 plot((1:length(b))/fs,b,'-b'); 0489 xlabel(['Time (' xticksi 's)']); 0490 if ixs>0 0491 title(spty{ixs}); 0492 end 0493 end 0494 if any(m=='S') && any(m1=='it') && numel(b)>0.1*fs % plot spectrogram 0495 if nfig 0496 figure(); 0497 end 0498 clf; 0499 if any(m=='L') 0500 sm='pJcwl'; 0501 else 0502 sm='pJcw'; 0503 end 0504 spgrambw(b,fs,sm); 0505 if ixs>0 0506 title(spty{ixs}); 0507 end 0508 end 0509 end