v_stdspectrum

PURPOSE ^

V_STDSPECTRUM Generate standard acoustic/speech spectra in s- or z-domain [B,A,SI,SN]=(S,M,F,N,ZI,BS,AS)

SYNOPSIS ^

function [b,a,si,sn]=v_stdspectrum(s,m,f,n,zi,bs,as)

DESCRIPTION ^

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]

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Generated by m2html © 2003