V_STDSPECTRUM Generate standard acoustic/speech spectra in s- or z-domain [B,A,SI,SN]=(S,M,F,N,ZI,BS,AS)
Usage: (1) [b,a]=v_stdspectrum(2,'z',fs) % create A-weighting z-domain filter for sample frequency fs
(2) [b,a]=v_stdspectrum(2,'zMLT',fs) % as above but also plot both s- and z-domain responses with log freq axis
(3) x=v_stdspectrum(11,'t',fs,n) % generate n samples of speech-shaped noise
(4) [b,a]=v_stdspectrum(0,'zEMLT',fs,[],[],bs,as) % convert s-domain filter bs(s)/as(s) to z-domain and plot approximation error
(5) for i=1:10; figure(i); v_stdspectrum(i,'z',3e4); end; v_tilefigs; % plot all the spectra for fs=30kHz
Inputs: s Spectrum type (either text or number - see below) or 0 to use bs/as
m mode: char 1 specifies output type (no combinations allowed),
f - frequency response (complex)
m - magnitude response
p - power spectrum
l - power per decade
d - decibel power spectrum
e - decibel power per decade
t - time waveform; the power equals the integral of the one-sided spectrum
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 in dB [default]
E - plot magnitude spectrum error in dB
F - magnitude spectrum error in dB with a -40 dB floor relative to peak
Q - plot phase spectrum
T - also plot target spectra [default]
A - plot zeros/poles
U - plot using a uniform frequency axis
L - plot using a log frequency axis [default]
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) if s=0
as denominator s-domain polynomial (or cell array containing polynomial factors if s=0
Outputs: b (1) numerator of the output spectrum (modes 's' or 'z')
(2) output waveform (mode 't')
(3) outptut spectrum (modes 'f', 'm', 'p' or 'd')
a (1) denonminator of the output spectrum (modes 's' or 'z')
(2) final state of the filter - use as the zi input of a future call (mode 't')
si spectrum type number (0 to 10)
sn spectrum name
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]
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].
11 LTASS-P50 : the long-term average speech spectrum taken from Table 1 of [4] in dB SPL @ 1m on-axis.
13 LTASS-1994 : the long-term average speech spectrum that is taken from Table 2 in [2] in dB SPL @ 1m on-axis.
14 EM3346-Gain : the response of a Knowles EM-3346 electret microphone in V/SPL from [12]
15 EM3346-Noise : the noise PSD of a Knowles EM-3346 electret microphone from [12]
Obsolete fits included for backward compatibility only:
5 X1-LTASS-P50 : (use 11 instead) the long-term average speech spectrum taken from Table 1 of [4].
6 X1-LTASS-1994 : (use 13 instead) the long-term average speech spectrum that is taken from Table 2 in [2]
12 X2-LTASS-1994 : (use 13 instead) the long-term average speech spectrum that is taken from Table 2 in [2]0001 function [b,a,si,sn]=v_stdspectrum(s,m,f,n,zi,bs,as) 0002 %V_STDSPECTRUM Generate standard acoustic/speech spectra in s- or z-domain [B,A,SI,SN]=(S,M,F,N,ZI,BS,AS) 0003 % 0004 % Usage: (1) [b,a]=v_stdspectrum(2,'z',fs) % create A-weighting z-domain filter for sample frequency fs 0005 % (2) [b,a]=v_stdspectrum(2,'zMLT',fs) % as above but also plot both s- and z-domain responses with log freq axis 0006 % (3) x=v_stdspectrum(11,'t',fs,n) % generate n samples of speech-shaped noise 0007 % (4) [b,a]=v_stdspectrum(0,'zEMLT',fs,[],[],bs,as) % convert s-domain filter bs(s)/as(s) to z-domain and plot approximation error 0008 % (5) for i=1:10; figure(i); v_stdspectrum(i,'z',3e4); end; v_tilefigs; % plot all the spectra for fs=30kHz 0009 % 0010 %Inputs: s Spectrum type (either text or number - see below) or 0 to use bs/as 0011 % m mode: char 1 specifies output type (no combinations allowed), 0012 % f - frequency response (complex) 0013 % m - magnitude response 0014 % p - power spectrum 0015 % l - power per decade 0016 % d - decibel power spectrum 0017 % e - decibel power per decade 0018 % t - time waveform; the power equals the integral of the one-sided spectrum 0019 % s - s-domain filter: b(s)/a(s) [default] 0020 % z - z-domain filter: b(z)/a(z) 0021 % i - sampled impulse response 0022 % plotting options 0023 % M - plot magnitude spectrum in dB [default] 0024 % E - plot magnitude spectrum error in dB 0025 % F - magnitude spectrum error in dB with a -40 dB floor relative to peak 0026 % Q - plot phase spectrum 0027 % T - also plot target spectra [default] 0028 % A - plot zeros/poles 0029 % U - plot using a uniform frequency axis 0030 % L - plot using a log frequency axis [default] 0031 % W - waveform 0032 % S - spectrogram 0033 % f sample frequency (modes 'z','i','t') or set of frequencies in Hz (modes 'f','m','p','d') 0034 % n number of output samples (mode 'i','t') 0035 % zi initial state of filter from a previous call (mode 't') 0036 % bs numerator s-domain polynomial (or cell array containing polynomial factors) if s=0 0037 % as denominator s-domain polynomial (or cell array containing polynomial factors if s=0 0038 % 0039 % Outputs: b (1) numerator of the output spectrum (modes 's' or 'z') 0040 % (2) output waveform (mode 't') 0041 % (3) outptut spectrum (modes 'f', 'm', 'p' or 'd') 0042 % a (1) denonminator of the output spectrum (modes 's' or 'z') 0043 % (2) final state of the filter - use as the zi input of a future call (mode 't') 0044 % si spectrum type number (0 to 10) 0045 % sn spectrum name 0046 % 0047 % Spectrum types (specify either as a number or case-insensitive text abbreviation): 0048 % 0 external : BS and AS arguments specify an s-domain filter 0049 % 1 White : white noise 0050 % 2 A-Weight : the formula for this is given in [3] and is based on 0051 % the equal-loudness curves of [9] 0052 % 3 B-Weight : this needs to be confirmed with ANSI S1.4-1981 standard or IEC 60651 0053 % 4 C-Weight : the formula for this is given in [3] 0054 % 7 SII-intinv : The inverse spectrum of the ear's internal masking noise; this is taken 0055 % from table 1 of [1]. It is inverted so that it is a bandpass rather than 0056 % bandstop characteristic. 0057 % 8 BS-468 : The weighting proposed for audio frequency noise measurement in [5] and [6]. 0058 % 9 USASI : Noise simulating long-term programme material spectrum from [7],[8]. 0059 % The level is such that the power is 0dB over an infinite bandwidth 0060 % 10 POTS : the D spectrum from [11]. 0061 % 11 LTASS-P50 : the long-term average speech spectrum taken from Table 1 of [4] in dB SPL @ 1m on-axis. 0062 % 13 LTASS-1994 : the long-term average speech spectrum that is taken from Table 2 in [2] in dB SPL @ 1m on-axis. 0063 % 14 EM3346-Gain : the response of a Knowles EM-3346 electret microphone in V/SPL from [12] 0064 % 15 EM3346-Noise : the noise PSD of a Knowles EM-3346 electret microphone from [12] 0065 % 0066 % Obsolete fits included for backward compatibility only: 0067 % 0068 % 5 X1-LTASS-P50 : (use 11 instead) the long-term average speech spectrum taken from Table 1 of [4]. 0069 % 6 X1-LTASS-1994 : (use 13 instead) the long-term average speech spectrum that is taken from Table 2 in [2] 0070 % 12 X2-LTASS-1994 : (use 13 instead) the long-term average speech spectrum that is taken from Table 2 in [2] 0071 0072 % References: 0073 % [1] Methods for the calculation of the speech intelligibility index. 0074 % ANSI Standard S3.5-1997 (R2007), American National Standards Institute, 1997. 0075 % [2] D. Byrne, H. Dillon, K. Tran, S. Arlinger, K. Wilbraham, R. Cox, B. Hayerman, 0076 % R. Hetu, J. Kei, C. Lui, J. Kiessling, M. N. Kotby, N. H. A. Nasser, 0077 % W. A. H. E. Kholy, Y. Nakanishi, H. Oyer, R. Powell, D. Stephens, R. Meredith, 0078 % T. Sirimanna, G. Tavartkiladze, G. I. Frolenkov, S. Westerman, and C. Ludvigsen. 0079 % An international comparison of long-term average speech spectra. 0080 % JASA, 96 (4): 2108-2120, Oct. 1994. 0081 % [3] CENELEC. Electroacoustics - sound level meters. Technical Report EN EN 61672-1:2003, 2003. 0082 % (also ANSI S1.42-2001) 0083 % [4] ITU-T. Artificial voices. Standard P.50, Sept. 1999. 0084 % [5] ITU-T. Measurement of weighted noise in sound-programme circuits. 0085 % Recommendation J.16, 1988. 0086 % [6] ITU-R. Measurement of audio-requency noise voltage level in sound broadcasting. 0087 % Recommendation BS.468., 1986. 0088 % [7] NRSC AM Reemphasis, Deemphasize, and Broadcast Audio Transmission Bandwidth Specifications, 0089 % EIA-549 Standard, Electronics Industries Association , July 1988. 0090 % [8] NRSC AM Reemphasis, Deemphasize, and Broadcast Audio Transmission Bandwidth Specifications, 0091 % NRSC-1-A Standard, Sept 2007, Online: http://www.nrscstandards.org/SG/NRSC-1-A.pdf 0092 % [9] H. Fletcher and W. A. Munson. Loudness, its definition, measurement and calculation. 0093 % J. Acoust Soc Amer, 5: 82-108, Oct. 1933. 0094 % [10] American National Standard Specification for Sound Level Meters. 0095 % ANSI S1.4-1983 (R2006)/ANSI S1.4a-1985 (R2006), American National Standards Institute 0096 % [11] IEEE standard equipment requirements and measurement techniques for analog transmission 0097 % parameters for telecommunications. Standard IEEE Std 743-1995, Dec. 1995. 0098 % [12] S.C.Thompson, J.L.LoPresti, E.M.Ring, H.G.Nepomuceno, J.J.Beard, W.J.Ballad, and E.V.Carlson. 0099 % Noise in miniature microphones. J. Acoust. Soc. Amer., 111 (2): 861�866, feb 2002. doi: 10.1121/1.1436072. 0100 0101 % Other candidates: (a) Z-weighting, (b) ISO226, (c) P.48 spectra 0102 % 0103 % Other standards: 0104 % IEEE743 has several weighting filters defined 0105 % ITU-T 0.41 Psophometer for use on telephone-type circuits 0106 % Bell System Technical Reference 41009 (C-message) 0107 % ISO 8041:2005 (E): Human Response to Vibration - Measuring Instrumentation 0108 % IEC 1260:1995, class 1 (also IEC 61260/ANSI S1.11-2004) Octave band and fractional octave band filters 0109 % IEC 651: Specification for Sound Level Meters 0110 % IRS P.48: sending and receiving characteristics defined by isolated points 0111 % mIRS P.830 modified IRS also defined by isolated points (see Annex D) available in G.191 0112 % G.191 software tools library contains IRS and mIRS implementations in FIR and IIR 0113 % ISO226 gives equal-loudness curves 0114 0115 % Copyright (C) Mike Brookes 2008-2018 0116 % Version: $Id: v_stdspectrum.m 10865 2018-09-21 17:22:45Z dmb $ 0117 % 0118 % VOICEBOX is a MATLAB toolbox for speech processing. 0119 % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 0120 % 0121 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0122 % This program is free software; you can redistribute it and/or modify 0123 % it under the terms of the GNU General Public License as published by 0124 % the Free Software Foundation; either version 2 of the License, or 0125 % (at your option) any later version. 0126 % 0127 % This program is distributed in the hope that it will be useful, 0128 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0129 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0130 % GNU General Public License for more details. 0131 % 0132 % You can obtain a copy of the GNU General Public License from 0133 % http://www.gnu.org/copyleft/gpl.html or by writing to 0134 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0135 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0136 0137 % Bugs/Suggestions: 0138 % * Could generate impulse response by appending a LP filter (elliptic) to 0139 % the s-domain transfer function and sampling the impulse response 0140 % * better calculation of impulse response length based on its total power 0141 % * ensure that the number of z-domain zeros at z=1 is correct. 0142 0143 persistent spty nspty spz ient fixz baz bazps spyax 0144 % spty contains the name of the spectrum 0145 % spyax contains the y axis label 0146 % spz contains a list of the poles and zeros 0147 % spz(1) = gain, spz(2) = number of zeros (excluding implicit conjugates), spz(3...) zeros followed by poles 0148 % only one of a complex conjugate pair is given. To print out b/a in the correct format use: 0149 % ra=roots(a);ra=ra(imag(ra)>=0);rb=roots(b);rb=rb(imag(rb)>=0);fprintf('\n[%0.14g %d',b(1)/a(1),length(rb)); 0150 % for i=1:length(rb), fprintf(' %s',v_sprintcpx(rb(i),'0.14gi')); end 0151 % for i=1:length(ra), fprintf(' %s',v_sprintcpx(ra(i),'0.14gi')); end, fprintf('];\n'); 0152 if isempty(spz) 0153 spty={'White';'A-Weight';'B-Weight';'C-Weight';'X1-LTASS-P50';'X1-LTASS-1994';'SII-IntInv';'BS-468';'USASI';'POTS';'LTASS-P50'; 0154 'X2-LTASS-1994';'LTASS-1994';'EM3346-Gain';'EM3346-Noise'}; 0155 nspty=size(spty,1); 0156 spyax={'';'';'';'';'dB-SPL/\surd{}Hz @ 25 mm';'';'SII-IntInv';'BS-468';'USASI';'POTS';'dB-SPL/\surd{}Hz @ 1m'; ''; 0157 'dB-SPL/\surd{}Hz @ 1m';'dB-V/SPL';'dB-V/\surd{}Hz'}; 0158 spz={[1 0]; 0159 [7390100803.6603 4 0 0 0 0 -129.42731565506 -129.42731565506 -676.40154023295 -4636.125126885 -76618.526016858 -76618.526016858]; 0160 [5986190155.0156 3 0 0 0 -129.42731565506 -129.42731565506 -995.88487118796 -76618.526016858 -76618.526016858]; 0161 [5912384617.784 2 0 0 -129.42731565506 -129.42731565506 -76618.526016858 -76618.526016858]; 0162 [1.1294790345421e+015 3 0 0 -34437.856184098 -721.94747118664+754.97798119504i -1721.704402273 -5234.2950286868 -10953.570773216+42789.342252749i]; 0163 [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]; 0164 [6.1311266394354e+018 2 0 0 -381.08293630892 -5920.974797779 -4701.76218192+24369.279310049i 10597.854874768+39258.915154617i]; 0165 [2.1034520039796e+024 1 0 -25903.701047817 -23615.535213635+36379.908937329i -62675.170058468 -18743.746690721+62460.156452506i]; 0166 [72.648989380657 1 0 -2*pi*[100 320]]; 0167 [7.8820088171767e+016 4 0 0 0 0 -452.681+1924.28i -2334+1702.73i -11264.2+8213.32i -4665.8+19828.7i]; 0168 [2.67e+013 3 0 0 -34437.856184098 -721.94747118664+754.97798119504i -1721.704402273 -5234.2950286868 -10953.570773216+42789.342252749i]; 0169 [1600000000 5 0 0 0 -22550.637578954 -11319.635610404+70239.177107659i -253.31327846696+672.10855509952i -1299.1885437118+2301.2064056419i -10646.952627978+68290.702816027i -147307.51763333 -628.3185307180 ]; 0170 [0.40294607181247 6 0 0 0 0 -103821.41495527 -6138.6617784378 -1387.7780129857+2694.482976041i -212.70001558505+701.9845877834i -489.69089908576 -241.48780111882]; 0171 [-2.309644968403e+23 2 0 -58249.19407694 -23164.744808779+191400.50379156i -34701.140307771+126665.52643837i -11829.456694949+32636.504600476i -33607.63134518 -991.35267796254]; 0172 [1.8926228601879e+14 4 -37153.159597416+720298.30892497i -6408.3219736031+187320.82306394i -37654.306739587+61707.365378187i -8055.8650530296 -22183.512456212+641448.91012113i -114814.21280417+529521.53918435i -93068.804147647+170656.9265202i -9742.2020045189+186595.22125487i -15084.001933885+34453.016304251i -898.60967839055]; 0173 }; 0174 bazps=cell(nspty,4); 0175 for i=1:nspty 0176 spzi=spz{i}; 0177 nsz=spzi(2); 0178 sz=spzi(3:3+nsz-1); 0179 sz=[sz conj(sz(imag(sz)~=0))]; 0180 sp=spzi(3+nsz:end); 0181 sp=[sp conj(sp(imag(sp)~=0))]; 0182 bazps{i,1}=spzi(1)*poly(sz); 0183 bazps{i,2}=poly(sp); 0184 bazps{i,3}=sz; 0185 bazps{i,4}=sp; 0186 end 0187 nz=15; % size of cache 0188 ient=0; % cache entry number 0189 fixz=repmat(-1,nz,2); 0190 baz=cell(nz,2); 0191 0192 end 0193 if nargin<2 || ~numel(m) 0194 m=' '; 0195 end 0196 m1=m(1); % output format 0197 if ~any(m1=='fmpldetszi') 0198 m1=char('s'+~nargout*('d'-'s')); % 's' normally, 'd' if no outputs 0199 end 0200 if nargin<3 0201 f=8192; % default frequency 0202 end 0203 0204 % determine the spectrum type 0205 0206 if ~numel(s) || s(1)==0 0207 si=0; 0208 sn=''; 0209 sb=1; 0210 sz=[]; % list of s-domain zeros 0211 if iscell(bs) 0212 for i=1:numel(bs) 0213 sb=conv(sb,bs{i}); 0214 sz=[sz roots(bs{i}).']; 0215 end 0216 elseif numel(bs)>0 0217 sb=bs; 0218 sz=roots(bs).'; 0219 end 0220 sa=1; 0221 sp=[]; % list of s-domain poles 0222 if iscell(as) 0223 for i=1:numel(as) 0224 sa=conv(sa,as{i}); 0225 sp=[sp roots(as{i}).']; 0226 end 0227 elseif numel(as)>0 0228 sa=as; 0229 sp=roots(as).'; 0230 end 0231 else 0232 if ischar(s) 0233 si=find(strcmpi(s,spty)); 0234 if isempty(si) 0235 error('undefined spectrum type: %s',s); 0236 end 0237 else 0238 si=s; 0239 end 0240 if si>nspty 0241 error('undefined spectrum type: %d',si); 0242 end 0243 sn=spty{si}; % name of spectrum 0244 % get s-domain function 0245 sb=bazps{si,1}; % sb/sa is transfer function 0246 sa=bazps{si,2}; 0247 sz=bazps{si,3}; % sz and sp are lists of zeros and poles 0248 sp=bazps{si,4}; 0249 end 0250 if (nargin<3 || ~numel(f)) && any(m1=='fmpd') % calcualte the frequency range 0251 apz=abs([sp sz]); 0252 apz(apz==0)=[]; % ignore zero frequency poles/zeros 0253 if ~numel(apz) 0254 apz=[100 5000]; 0255 elseif length(apz)==1 0256 apz=[apz/10 apz*10]; 0257 end 0258 f=logspace(log10(min(apz)*0.5/pi)-0.5,log10(max(apz)*0.5/pi)+0.5,200); 0259 end 0260 if any(m1=='fmpdle') 0261 h=polyval(sb,2i*pi*f)./polyval(sa,2i*pi*f); 0262 end 0263 fs=f; % save sampling frequency 0264 if any(m1=='izt') % we need a z-domain filter 0265 if si==1 % treat white noise specially 0266 bz=1; 0267 az=1; 0268 else 0269 jent=find(fixz(:,2)==f*nspty+si,1); % see if it is in the cache already 0270 if numel(jent) 0271 ient=ient+1; % cache entry number 0272 fixz(jent,1)=ient; % update index to show it is recently used 0273 bz=baz{jent,1}; 0274 az=baz{jent,2}; 0275 else 0276 % we use an iterative method to find the best digital filter 0277 % we initialize the phases with either bilinear or impulse invariance 0278 % only using the impulse invariance if it is good (max error < 10dB) 0279 % we then iterate invfreqz using the s-domain magnitudes and 0280 % the phases of the best fit so far. 0281 % we use log-spaced frequencies at low frequencies and linear at high 0282 % we then search for various numbers of poles and zeros 0283 warning off all % avoid lots of ill-conditioning error messages 0284 nflin=100; % number of frequency samples in linear region (high freq) 0285 alp=1.15; % freq ratio increment in low freq region 0286 f0=10*2*pi/f; % minimum interesting frequency (10 Hz in radians) 0287 fx=pi/nflin/(alp-1); % boundary between log and linear portions 0288 if fx<=f0 || f0>=pi 0289 fif=linspace(0,pi,nflin); 0290 elseif fx>pi 0291 fif=[0 logspace(log10(f0),log10(pi),ceil(log10(pi/f0)/log10(alp)))]; 0292 else 0293 nlin=ceil((pi-fx)*nflin/pi); 0294 fif=[0 logspace(log10(f0),log10(fx),ceil(log10(fx/f0)/log10(alp))) linspace(fx+(pi-fx)/nlin,pi,nlin-1)]; 0295 end 0296 h0=abs(polyval(sb,1i*fif*f)./polyval(sa,1i*fif*f)); % target magnitude spectrum 0297 h0tol=max(h0)*1e-4; % absolute gain tolerance 0298 % initialize with impulse invariance 0299 [bz,az]=impinvar(sb,sa,f); 0300 hj=freqz(bz,az,fif); 0301 maxj=max(abs(db(abs(hj)+h0tol)-db(h0+h0tol))); 0302 % or else with bilinear 0303 [ifb,ifa]=bilinear(sb,sa,f); 0304 hn=freqz(ifb,ifa,fif); 0305 maxi=max(abs(db(abs(hn+h0tol))-db(h0+h0tol))); 0306 if maxi<maxj || maxj>10 % accept bilinear if it is better or if imp inv is bad 0307 maxj=maxi; 0308 bz=ifb; 0309 az=ifa; 0310 hj=hn; 0311 end 0312 pat0=sb(end)==0; % we have a zero at DC 0313 if pat0 0314 fif(1)=[]; % eliminate DC as a probe frequency 0315 hz1=1-exp(-1i*fif); % response of zero at z=1 0316 h0=h0(2:end)./abs(hz1); % remove a zero at z=1 from the target 0317 hj=hj(2:end)./hz1; % remove a zero at z=1 from the initial phase 0318 end 0319 upd=0; 0320 for mm=1:length(sa) % maximum number of poles 0321 for nn=1:mm % number of zeros is always less than number of poles 0322 hn=hj; 0323 j=0; 0324 for i=1:30 % iterate up to 30 times (usually less) 0325 h=h0.*exp(1i*angle(hn)); 0326 [ifb,ifa]=invfreqz(h,fif,nn,mm,[],10); 0327 hn=freqz(ifb,ifa,fif); 0328 maxi=max(abs(db(abs(hn+h0tol))-db(h0+h0tol))); 0329 if maxi<maxj 0330 maxj=maxi; 0331 bz=ifb; 0332 az=ifa; 0333 hj=hn; 0334 j=i; 0335 upd=1; 0336 end 0337 if i>j+5 % quit if no improvement in last five iterations 0338 break 0339 end 0340 end 0341 end 0342 end 0343 if upd 0344 bz=conv(bz,[1 -1]); % restore the zero at z=0 0345 end 0346 warning on all 0347 if si>0 0348 ient=ient+1; % cache entry number 0349 [jdum,jent]=min(fixz(:,1)); % find least recently used cache entry 0350 fixz(jent,1)=ient; % flag as recently used 0351 fixz(jent,2)=f*nspty+si; % save frequency/spectrum code 0352 baz{jent,1}=bz; 0353 baz{jent,2}=az; 0354 end 0355 end 0356 end 0357 end 0358 switch m1 0359 case 'z' 0360 b=bz; 0361 a=az; 0362 case 't' 0363 if nargin<5 || ~numel(zi) 0364 [b,a]=v_randfilt(bz,az,n); 0365 else 0366 [b,a]=v_randfilt(bz,az,n,zi); 0367 end 0368 b=b*10*log10(fs/2); % scale to compensate for bandwidth 0369 case 'i' 0370 if nargin<5 || ~numel(zi) 0371 if nargin<4 || ~numel(n) % determine n to include 1 - 1e-8 of the energy 0372 n=ceil(-fs*log(1e4)/max(real(sp))); 0373 end 0374 [b,a]=filter(bz,az,[1; zeros(n-1,1)]); 0375 else 0376 [b,a]=filter(bz,az,zeros(n,1),zi); 0377 end 0378 case 'm' 0379 b = abs(h); 0380 case 'f' 0381 b = h; 0382 case 'd' 0383 b = db(abs(h)); 0384 case 'e' 0385 b=db(abs(h).*f*log(10)); % convert to power per decade in dB 0386 case 'l' 0387 b=h.*conj(h).*f*log(10); % convert to power per decade 0388 case 'p' 0389 b=h.*conj(h); 0390 case 's' 0391 b=sb; 0392 a=sa; 0393 otherwise 0394 error('Output format %s not implemented',m1); 0395 end 0396 0397 % plot data 0398 if ~nargout || ~strcmp(m,lower(m)) 0399 if strcmp(m,lower(m)) % if no upper case letters 0400 m='MLT'; 0401 end 0402 if ~any(m=='Q') && ~any(m=='A') && ~any(m=='E') && ~any(m=='F') && (m1~='t' || ~any(m=='W') && ~any(m=='S')) 0403 m(end+1)='M'; % default plot type 0404 end 0405 nfig=0; 0406 paz=any(m1=='itz'); % plot discrete time result 0407 pas=~paz || any(m=='T'); % plot continuous time result 0408 if any(m=='M') || any(m=='Q') || any(m=='E') || any(m=='E') % draw a frequency response plot 0409 clf; 0410 nfig=1; 0411 pam=any(m=='M'); % magnitude response 0412 pae=(any(m=='E') || any(m=='F')) && paz; % magnitude response error 0413 paq=any(m=='Q'); % phase response 0414 pat=any(m=='T') && paz; % include target spectrum 0415 pal=any(m=='L') || (~paz && ~any(m=='U')); % log frequency axis 0416 if any(m1=='itz') 0417 fs=f; % save the sample frequency 0418 apz=abs([sp sz]); % and determine the frequency range 0419 apz(apz==0)=[]; % ignore zero frequency poles/zeros 0420 if ~numel(apz) 0421 apz=[100 5000]; 0422 elseif length(apz)==1 0423 apz=[apz/10 apz*10]; 0424 end 0425 if pal 0426 f=logspace(log10(min([fs/1000 apz*0.05/pi])),log10(fs/2),200); 0427 else 0428 f=linspace(min([fs/1000 apz*0.05/pi]),fs/2,200); 0429 end 0430 end 0431 hs=freqs(sb,sa,2*pi*f); % s domain response 0432 if paz 0433 hz=freqz(bz,az,f,fs); % z domain response 0434 end 0435 axh=[]; 0436 nax=pam+pae+paq; % number of axis sets 0437 titex=''; 0438 if pam % plot a magnitude response 0439 nf=length(f); 0440 df=0.5*(f([2:nf nf])-f([1 1:nf-1]))'; 0441 if nax>1 0442 subplot(nax,1,1); 0443 end 0444 if paz 0445 if pas % plot s domain and z domain magnitude plots 0446 plot(f,db(abs(hs)),'--r',f,db(abs(hz)),'-b') 0447 ymax=max(db(abs([hs hz])))+1; 0448 ymin=min(db(abs([hs hz])))-1; 0449 titex=' ( : = s, - = z)'; 0450 pwrint=sprintf('\\int=%.1f, %.1f dB',db(abs(hs).^2*df)/2,db(abs(hz).^2*df)/2); 0451 % legend('s-domain','z-domain','location','best'); 0452 else % plot z domain magnitude plot only 0453 plot(f,db(abs(hz)),'-b') 0454 ymax=max(db(abs(hz)))+1; 0455 ymin=min(db(abs(hz)))-1; 0456 pwrint=sprintf('\\int=%.1f dB',db(abs(hz).^2*df)/2); 0457 end 0458 else % plot s domain magnitude plot only 0459 plot(f,db(abs(hs)),'-b') 0460 ymax=max(db(abs(hs)))+1; 0461 ymin=min(db(abs(hs)))-1; 0462 pwrint=sprintf('\\int=%.1f dB',db(abs(hs).^2*df)/2); 0463 end 0464 if pal 0465 set(gca,'Xscale','log'); 0466 end 0467 set(gca,'YLim',[max(ymin,ymax-60) ymax]); 0468 v_axisenlarge([-1 1.05]); 0469 v_texthvc(0.02,0.98,pwrint,'LTk'); 0470 xlabel(['Frequency (' v_xticksi 'Hz)']); 0471 if si>0 && ~isempty(spyax{si}) 0472 ylabel(spyax{si}); 0473 else 0474 ylabel('Gain (dB)'); 0475 end 0476 0477 if si>0 0478 title(sprintf('Type %d: %s%s',si,spty{si},titex)); 0479 end 0480 axh(end+1)=gca; 0481 end 0482 if pae % plot magnitude response error 0483 if nax>1 0484 subplot(nax,1,1+pam); 0485 end 0486 if any(m=='F') 0487 dbflr=40; 0488 dbfl=max(abs(hs))*10^(-dbflr/20); % make a floor 40 dB below peak 0489 dbflt=sprintf(' (floor@-%.0fdB)',dbflr); 0490 else 0491 dbfl=0; 0492 dbflt=''; 0493 end 0494 dberr=db(abs(hz)+dbfl)-db(abs(hs)+dbfl); 0495 plot([f(1) f(end)],[0 0],':k',f,dberr,'-b') 0496 if pal 0497 set(gca,'Xscale','log'); 0498 end 0499 v_axisenlarge([-1 -1.05]); 0500 % set(gca,'XLim',[min(f) max(f)]); 0501 xlabel(['Frequency (' v_xticksi 'Hz)']); 0502 ylabel('Gain Error (dB)'); 0503 v_texthvc(0.95,0.95,sprintf('Err < %.1f dB%s',max(abs(dberr)),dbflt),'RTk'); 0504 if si>0 && ~pam 0505 title(sprintf('Type %d: %s%s',si,spty{si},titex)); 0506 end 0507 axh(end+1)=gca; 0508 end 0509 if paq 0510 if nax>1 0511 subplot(nax,1,nax); 0512 end 0513 if paz 0514 if pas 0515 plot(f,angle(hs),'--r',f,angle(hz),'-b') 0516 else 0517 plot(f,angle(hz),'-b') 0518 end 0519 else 0520 plot(f,angle(hs),'-b') 0521 end 0522 if pal 0523 set(gca,'Xscale','log'); 0524 end 0525 v_axisenlarge([-1 -1.05]); 0526 % set(gca,'XLim',[min(f) max(f)]); 0527 xlabel(['Frequency (' v_xticksi 'Hz)']); 0528 ylabel('Phase (rad)'); 0529 if si>0 && nax==1 0530 title(spty{si}); 0531 end 0532 axh(end+1)=gca; 0533 end 0534 if nax>1 0535 linkaxes(axh,'x'); 0536 end 0537 end 0538 if any(m=='A') % plot complex plane 0539 if nfig 0540 figure(); 0541 end 0542 clf; 0543 nfig=1; 0544 if pas 0545 if paz 0546 subplot(121); 0547 end 0548 plot(real(sp),imag(sp),'xb',real(sz),imag(sz),'ob'); 0549 axis equal; 0550 xlim=get(gca,'xlim'); 0551 xlim(1)=min(xlim(1),-1000); 0552 xlim(2)=max(xlim(2),1000); 0553 ylim=get(gca,'ylim'); 0554 ylim(1)=min(ylim(1),-1000); 0555 ylim(2)=max(ylim(2),1000); 0556 axis([xlim ylim]); 0557 hold on 0558 plot(xlim,[0 0],':r',[0 0],ylim,':r'); 0559 hold off 0560 title('s-plane'); 0561 end 0562 if paz 0563 if pas 0564 subplot(122); 0565 end 0566 axl=max(abs([1.1; az(:);bz(:)])); 0567 t=linspace(0,2*pi); 0568 rtzb=roots(bz); 0569 rtza=roots(az); 0570 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'); 0571 axis equal; 0572 axis([-1 1 -1 1]*axl); 0573 title('z-plane'); 0574 end 0575 end 0576 if any(m=='W') && any(m1=='it') % plot waveform 0577 if nfig 0578 figure(); 0579 end 0580 clf; 0581 nfig=1; 0582 plot((1:length(b))/fs,b,'-b'); 0583 xlabel(['Time (' v_xticksi 's)']); 0584 if si>0 0585 title(spty{si}); 0586 end 0587 end 0588 if any(m=='S') && any(m1=='it') && numel(b)>0.1*fs % plot spectrogram 0589 if nfig 0590 figure(); 0591 end 0592 clf; 0593 if any(m=='L') 0594 sm='pJcwl'; 0595 else 0596 sm='pJcw'; 0597 end 0598 v_spgrambw(b,fs,sm); 0599 if si>0 0600 title(spty{si}); 0601 end 0602 end 0603 end