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 % * have option to return the filter as biquads for better stability at high sample rates 0143 0144 persistent spty nspty spz ient fixz baz bazps spyax 0145 % spty contains the name of the spectrum 0146 % spyax contains the y axis label 0147 % spz contains a list of the poles and zeros 0148 % spz(1) = gain, spz(2) = number of zeros (excluding implicit conjugates), spz(3...) zeros followed by poles 0149 % only one of a complex conjugate pair is given. To print out b/a in the correct format use: 0150 % 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)); 0151 % for i=1:length(rb), fprintf(' %s',v_sprintcpx(rb(i),'0.14gi')); end 0152 % for i=1:length(ra), fprintf(' %s',v_sprintcpx(ra(i),'0.14gi')); end, fprintf('];\n'); 0153 if isempty(spz) 0154 spty={'White';'A-Weight';'B-Weight';'C-Weight';'X1-LTASS-P50';'X1-LTASS-1994';'SII-IntInv';'BS-468';'USASI';'POTS';'LTASS-P50'; 0155 'X2-LTASS-1994';'LTASS-1994';'EM3346-Gain';'EM3346-Noise'}; 0156 nspty=size(spty,1); 0157 spyax={'';'';'';'';'dB-SPL/\surd{}Hz @ 25 mm';'';'SII-IntInv';'BS-468';'USASI';'POTS';'dB-SPL/\surd{}Hz @ 1m'; ''; 0158 'dB-SPL/\surd{}Hz @ 1m';'dB-V/SPL';'dB-V/\surd{}Hz'}; 0159 spz={[1 0]; 0160 [7390100803.6603 4 0 0 0 0 -129.42731565506 -129.42731565506 -676.40154023295 -4636.125126885 -76618.526016858 -76618.526016858]; 0161 [5986190155.0156 3 0 0 0 -129.42731565506 -129.42731565506 -995.88487118796 -76618.526016858 -76618.526016858]; 0162 [5912384617.784 2 0 0 -129.42731565506 -129.42731565506 -76618.526016858 -76618.526016858]; 0163 [1.1294790345421e+015 3 0 0 -34437.856184098 -721.94747118664+754.97798119504i -1721.704402273 -5234.2950286868 -10953.570773216+42789.342252749i]; 0164 [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]; 0165 [6.1311266394354e+018 2 0 0 -381.08293630892 -5920.974797779 -4701.76218192+24369.279310049i 10597.854874768+39258.915154617i]; 0166 [2.1034520039796e+024 1 0 -25903.701047817 -23615.535213635+36379.908937329i -62675.170058468 -18743.746690721+62460.156452506i]; 0167 [72.648989380657 1 0 -2*pi*[100 320]]; 0168 [7.8820088171767e+016 4 0 0 0 0 -452.681+1924.28i -2334+1702.73i -11264.2+8213.32i -4665.8+19828.7i]; 0169 [2.67e+013 3 0 0 -34437.856184098 -721.94747118664+754.97798119504i -1721.704402273 -5234.2950286868 -10953.570773216+42789.342252749i]; 0170 [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 ]; 0171 [0.40294607181247 6 0 0 0 0 -103821.41495527 -6138.6617784378 -1387.7780129857+2694.482976041i -212.70001558505+701.9845877834i -489.69089908576 -241.48780111882]; 0172 [-2.309644968403e+23 2 0 -58249.19407694 -23164.744808779+191400.50379156i -34701.140307771+126665.52643837i -11829.456694949+32636.504600476i -33607.63134518 -991.35267796254]; 0173 [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]; 0174 }; 0175 bazps=cell(nspty,4); 0176 for i=1:nspty 0177 spzi=spz{i}; 0178 nsz=spzi(2); 0179 sz=spzi(3:3+nsz-1); 0180 sz=[sz conj(sz(imag(sz)~=0))]; 0181 sp=spzi(3+nsz:end); 0182 sp=[sp conj(sp(imag(sp)~=0))]; 0183 bazps{i,1}=spzi(1)*poly(sz); 0184 bazps{i,2}=poly(sp); 0185 bazps{i,3}=sz; 0186 bazps{i,4}=sp; 0187 end 0188 nz=15; % size of cache 0189 ient=0; % cache entry number 0190 fixz=repmat(-1,nz,2); 0191 baz=cell(nz,2); 0192 0193 end 0194 if nargin<2 || ~numel(m) 0195 m=' '; 0196 end 0197 m1=m(1); % output format 0198 if ~any(m1=='fmpldetszi') 0199 m1=char('s'+~nargout*('d'-'s')); % 's' normally, 'd' if no outputs 0200 end 0201 if nargin<3 0202 f=8192; % default frequency 0203 end 0204 0205 % determine the spectrum type 0206 0207 if ~numel(s) || s(1)==0 0208 si=0; 0209 sn=''; 0210 sb=1; 0211 sz=[]; % list of s-domain zeros 0212 if iscell(bs) 0213 for i=1:numel(bs) 0214 sb=conv(sb,bs{i}); 0215 sz=[sz roots(bs{i}).']; 0216 end 0217 elseif numel(bs)>0 0218 sb=bs; 0219 sz=roots(bs).'; 0220 end 0221 sa=1; 0222 sp=[]; % list of s-domain poles 0223 if iscell(as) 0224 for i=1:numel(as) 0225 sa=conv(sa,as{i}); 0226 sp=[sp roots(as{i}).']; 0227 end 0228 elseif numel(as)>0 0229 sa=as; 0230 sp=roots(as).'; 0231 end 0232 else 0233 if ischar(s) 0234 si=find(strcmpi(s,spty)); 0235 if isempty(si) 0236 error('undefined spectrum type: %s',s); 0237 end 0238 else 0239 si=s; 0240 end 0241 if si>nspty 0242 error('undefined spectrum type: %d',si); 0243 end 0244 sn=spty{si}; % name of spectrum 0245 % get s-domain function 0246 sb=bazps{si,1}; % sb/sa is transfer function 0247 sa=bazps{si,2}; 0248 sz=bazps{si,3}; % sz and sp are lists of zeros and poles 0249 sp=bazps{si,4}; 0250 end 0251 if (nargin<3 || ~numel(f)) && any(m1=='fmpd') % calcualte the frequency range 0252 apz=abs([sp sz]); 0253 apz(apz==0)=[]; % ignore zero frequency poles/zeros 0254 if ~numel(apz) 0255 apz=[100 5000]; 0256 elseif length(apz)==1 0257 apz=[apz/10 apz*10]; 0258 end 0259 f=logspace(log10(min(apz)*0.5/pi)-0.5,log10(max(apz)*0.5/pi)+0.5,200); 0260 end 0261 if any(m1=='fmpdle') 0262 h=polyval(sb,2i*pi*f)./polyval(sa,2i*pi*f); 0263 end 0264 fs=f; % save sampling frequency 0265 if any(m1=='izt') % we need a z-domain filter 0266 if si==1 % treat white noise specially 0267 bz=1; 0268 az=1; 0269 else 0270 jent=find(fixz(:,2)==f*nspty+si,1); % see if it is in the cache already 0271 if numel(jent) 0272 ient=ient+1; % cache entry number 0273 fixz(jent,1)=ient; % update index to show it is recently used 0274 bz=baz{jent,1}; 0275 az=baz{jent,2}; 0276 else 0277 % we use an iterative method to find the best digital filter 0278 % we initialize the phases with either bilinear or impulse invariance 0279 % only using the impulse invariance if it is good (max error < 10dB) 0280 % we then iterate invfreqz using the s-domain magnitudes and 0281 % the phases of the best fit so far. 0282 % we use log-spaced frequencies at low frequencies and linear at high 0283 % we then search for various numbers of poles and zeros 0284 warning off all % avoid lots of ill-conditioning error messages 0285 nflin=100; % number of frequency samples in linear region (high freq) 0286 alp=1.15; % freq ratio increment in low freq region 0287 f0=10*2*pi/f; % minimum interesting frequency (10 Hz in radians) 0288 fx=pi/nflin/(alp-1); % boundary between log and linear portions 0289 if fx<=f0 || f0>=pi 0290 fif=linspace(0,pi,nflin); 0291 elseif fx>pi 0292 fif=[0 logspace(log10(f0),log10(pi),ceil(log10(pi/f0)/log10(alp)))]; 0293 else 0294 nlin=ceil((pi-fx)*nflin/pi); 0295 fif=[0 logspace(log10(f0),log10(fx),ceil(log10(fx/f0)/log10(alp))) linspace(fx+(pi-fx)/nlin,pi,nlin-1)]; 0296 end 0297 h0=abs(polyval(sb,1i*fif*f)./polyval(sa,1i*fif*f)); % target magnitude spectrum 0298 h0tol=max(h0)*1e-4; % absolute gain tolerance 0299 % initialize with impulse invariance 0300 [bz,az]=impinvar(sb,sa,f); 0301 hj=freqz(bz,az,fif); 0302 maxj=max(abs(db(abs(hj)+h0tol)-db(h0+h0tol))); 0303 % or else with bilinear 0304 [ifb,ifa]=bilinear(sb,sa,f); 0305 hn=freqz(ifb,ifa,fif); 0306 maxi=max(abs(db(abs(hn+h0tol))-db(h0+h0tol))); 0307 if maxi<maxj || maxj>10 % accept bilinear if it is better or if imp inv is bad 0308 maxj=maxi; 0309 bz=ifb; 0310 az=ifa; 0311 hj=hn; 0312 end 0313 pat0=sb(end)==0; % we have a zero at DC 0314 if pat0 0315 fif(1)=[]; % eliminate DC as a probe frequency 0316 hz1=1-exp(-1i*fif); % response of zero at z=1 0317 h0=h0(2:end)./abs(hz1); % remove a zero at z=1 from the target 0318 hj=hj(2:end)./hz1; % remove a zero at z=1 from the initial phase 0319 end 0320 upd=0; 0321 for mm=1:length(sa) % maximum number of poles 0322 for nn=1:mm % number of zeros is always less than number of poles 0323 hn=hj; 0324 j=0; 0325 for i=1:30 % iterate up to 30 times (usually less) 0326 h=h0.*exp(1i*angle(hn)); 0327 [ifb,ifa]=invfreqz(h,fif,nn,mm,[],10); 0328 hn=freqz(ifb,ifa,fif); 0329 maxi=max(abs(db(abs(hn+h0tol))-db(h0+h0tol))); 0330 if maxi<maxj 0331 maxj=maxi; 0332 bz=ifb; 0333 az=ifa; 0334 hj=hn; 0335 j=i; 0336 upd=1; 0337 end 0338 if i>j+5 % quit if no improvement in last five iterations 0339 break 0340 end 0341 end 0342 end 0343 end 0344 if upd 0345 bz=conv(bz,[1 -1]); % restore the zero at z=0 0346 end 0347 warning on all 0348 if si>0 0349 ient=ient+1; % cache entry number 0350 [jdum,jent]=min(fixz(:,1)); % find least recently used cache entry 0351 fixz(jent,1)=ient; % flag as recently used 0352 fixz(jent,2)=f*nspty+si; % save frequency/spectrum code 0353 baz{jent,1}=bz; 0354 baz{jent,2}=az; 0355 end 0356 end 0357 end 0358 end 0359 switch m1 0360 case 'z' 0361 b=bz; 0362 a=az; 0363 case 't' 0364 if nargin<5 || ~numel(zi) 0365 [b,a]=v_randfilt(bz,az,n); 0366 else 0367 [b,a]=v_randfilt(bz,az,n,zi); 0368 end 0369 b=b*10*log10(fs/2); % scale to compensate for bandwidth 0370 case 'i' 0371 if nargin<5 || ~numel(zi) 0372 if nargin<4 || ~numel(n) % determine n to include 1 - 1e-8 of the energy 0373 n=ceil(-fs*log(1e4)/max(real(sp))); 0374 end 0375 [b,a]=filter(bz,az,[1; zeros(n-1,1)]); 0376 else 0377 [b,a]=filter(bz,az,zeros(n,1),zi); 0378 end 0379 case 'm' 0380 b = abs(h); 0381 case 'f' 0382 b = h; 0383 case 'd' 0384 b = db(abs(h)); 0385 case 'e' 0386 b=db(abs(h).*f*log(10)); % convert to power per decade in dB 0387 case 'l' 0388 b=h.*conj(h).*f*log(10); % convert to power per decade 0389 case 'p' 0390 b=h.*conj(h); 0391 case 's' 0392 b=sb; 0393 a=sa; 0394 otherwise 0395 error('Output format %s not implemented',m1); 0396 end 0397 0398 % plot data 0399 if ~nargout || ~strcmp(m,lower(m)) 0400 if strcmp(m,lower(m)) % if no upper case letters 0401 m='MLT'; 0402 end 0403 if ~any(m=='Q') && ~any(m=='A') && ~any(m=='E') && ~any(m=='F') && (m1~='t' || ~any(m=='W') && ~any(m=='S')) 0404 m(end+1)='M'; % default plot type 0405 end 0406 nfig=0; 0407 paz=any(m1=='itz'); % plot discrete time result 0408 pas=~paz || any(m=='T'); % plot continuous time result 0409 if any(m=='M') || any(m=='Q') || any(m=='E') || any(m=='E') % draw a frequency response plot 0410 clf; 0411 nfig=1; 0412 pam=any(m=='M'); % magnitude response 0413 pae=(any(m=='E') || any(m=='F')) && paz; % magnitude response error 0414 paq=any(m=='Q'); % phase response 0415 pat=any(m=='T') && paz; % include target spectrum 0416 pal=any(m=='L') || (~paz && ~any(m=='U')); % log frequency axis 0417 if any(m1=='itz') 0418 fs=f; % save the sample frequency 0419 apz=abs([sp sz]); % and determine the frequency range 0420 apz(apz==0)=[]; % ignore zero frequency poles/zeros 0421 if ~numel(apz) 0422 apz=[100 5000]; 0423 elseif length(apz)==1 0424 apz=[apz/10 apz*10]; 0425 end 0426 if pal 0427 f=logspace(log10(min([fs/1000 apz*0.05/pi])),log10(fs/2),200); 0428 else 0429 f=linspace(min([fs/1000 apz*0.05/pi]),fs/2,200); 0430 end 0431 end 0432 hs=freqs(sb,sa,2*pi*f); % s domain response 0433 if paz 0434 hz=freqz(bz,az,f,fs); % z domain response 0435 end 0436 axh=[]; 0437 nax=pam+pae+paq; % number of axis sets 0438 titex=''; 0439 if pam % plot a magnitude response 0440 nf=length(f); 0441 df=0.5*(f([2:nf nf])-f([1 1:nf-1]))'; 0442 if nax>1 0443 subplot(nax,1,1); 0444 end 0445 if paz 0446 if pas % plot both s domain and z domain magnitude plots 0447 plot(f,db(abs(hs)),'--r',f,db(abs(hz)),'-b') 0448 ymax=max(db(abs([hs hz])))+1; 0449 ymin=min(db(abs([hs hz])))-1; 0450 titex=' ( : = s, - = z)'; 0451 pwrint=sprintf('\\int df = %.1f, %.1f dB',db(abs(hs).^2*df)/2,db(abs(hz).^2*df)/2); 0452 % legend('s-domain','z-domain','location','best'); 0453 else % plot z domain magnitude plot only 0454 plot(f,db(abs(hz)),'-b') 0455 ymax=max(db(abs(hz)))+1; 0456 ymin=min(db(abs(hz)))-1; 0457 pwrint=sprintf('\\int df = %.1f dB',db(abs(hz).^2*df)/2); 0458 end 0459 else % plot s domain magnitude plot only 0460 plot(f,db(abs(hs)),'-b') 0461 ymax=max(db(abs(hs)))+1; 0462 ymin=min(db(abs(hs)))-1; 0463 pwrint=sprintf('\\int df = %.1f dB',db(abs(hs).^2*df)/2); 0464 end 0465 if pal 0466 set(gca,'Xscale','log'); 0467 end 0468 set(gca,'YLim',[max(ymin,ymax-60) ymax]); 0469 v_axisenlarge([-1 1.05]); 0470 v_texthvc(0.02,0.98,pwrint,'LTk'); 0471 xlabel(['Frequency (' v_xticksi 'Hz)']); 0472 if si>0 && ~isempty(spyax{si}) 0473 ylabel(spyax{si}); 0474 else 0475 ylabel('Gain (dB)'); 0476 end 0477 0478 if si>0 0479 title(sprintf('Type %d: %s%s',si,spty{si},titex)); 0480 end 0481 axh(end+1)=gca; 0482 end 0483 if pae % plot magnitude response error 0484 if nax>1 0485 subplot(nax,1,1+pam); 0486 end 0487 if any(m=='F') 0488 dbflr=40; 0489 dbfl=max(abs(hs))*10^(-dbflr/20); % make a floor 40 dB below peak 0490 dbflt=sprintf(' (floor@-%.0fdB)',dbflr); 0491 else 0492 dbfl=0; 0493 dbflt=''; 0494 end 0495 dberr=db(abs(hz)+dbfl)-db(abs(hs)+dbfl); 0496 plot([f(1) f(end)],[0 0],':k',f,dberr,'-b') 0497 if pal 0498 set(gca,'Xscale','log'); 0499 end 0500 v_axisenlarge([-1 -1.05]); 0501 % set(gca,'XLim',[min(f) max(f)]); 0502 xlabel(['Frequency (' v_xticksi 'Hz)']); 0503 ylabel('Gain Error (dB)'); 0504 v_texthvc(0.95,0.95,sprintf('Err < %.1f dB%s',max(abs(dberr)),dbflt),'RTk'); 0505 if si>0 && ~pam 0506 title(sprintf('Type %d: %s%s',si,spty{si},titex)); 0507 end 0508 axh(end+1)=gca; 0509 end 0510 if paq 0511 if nax>1 0512 subplot(nax,1,nax); 0513 end 0514 if paz 0515 if pas 0516 plot(f,angle(hs),'--r',f,angle(hz),'-b') 0517 else 0518 plot(f,angle(hz),'-b') 0519 end 0520 else 0521 plot(f,angle(hs),'-b') 0522 end 0523 if pal 0524 set(gca,'Xscale','log'); 0525 end 0526 v_axisenlarge([-1 -1.05]); 0527 % set(gca,'XLim',[min(f) max(f)]); 0528 xlabel(['Frequency (' v_xticksi 'Hz)']); 0529 ylabel('Phase (rad)'); 0530 if si>0 && nax==1 0531 title(spty{si}); 0532 end 0533 axh(end+1)=gca; 0534 end 0535 if nax>1 0536 linkaxes(axh,'x'); 0537 end 0538 end 0539 if any(m=='A') % plot complex plane 0540 if nfig 0541 figure(); 0542 end 0543 clf; 0544 nfig=1; 0545 if pas 0546 if paz 0547 subplot(121); 0548 end 0549 plot(real(sp),imag(sp),'xb',real(sz),imag(sz),'ob'); 0550 axis equal; 0551 xlim=get(gca,'xlim'); 0552 xlim(1)=min(xlim(1),-1000); 0553 xlim(2)=max(xlim(2),1000); 0554 ylim=get(gca,'ylim'); 0555 ylim(1)=min(ylim(1),-1000); 0556 ylim(2)=max(ylim(2),1000); 0557 axis([xlim ylim]); 0558 hold on 0559 plot(xlim,[0 0],':r',[0 0],ylim,':r'); 0560 hold off 0561 title('s-plane'); 0562 end 0563 if paz 0564 if pas 0565 subplot(122); 0566 end 0567 axl=max(abs([1.1; az(:);bz(:)])); 0568 t=linspace(0,2*pi); 0569 rtzb=roots(bz); 0570 rtza=roots(az); 0571 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'); 0572 axis equal; 0573 axis([-1 1 -1 1]*axl); 0574 title('z-plane'); 0575 end 0576 end 0577 if any(m=='W') && any(m1=='it') % plot waveform 0578 if nfig 0579 figure(); 0580 end 0581 clf; 0582 nfig=1; 0583 plot((1:length(b))/fs,b,'-b'); 0584 xlabel(['Time (' v_xticksi 's)']); 0585 if si>0 0586 title(spty{si}); 0587 end 0588 end 0589 if any(m=='S') && any(m1=='it') && numel(b)>0.1*fs % plot spectrogram 0590 if nfig 0591 figure(); 0592 end 0593 clf; 0594 if any(m=='L') 0595 sm='pJcwl'; 0596 else 0597 sm='pJcw'; 0598 end 0599 v_spgrambw(b,fs,sm); 0600 if si>0 0601 title(spty{si}); 0602 end 0603 end 0604 end