# 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:

- v_axisenlarge V_AXISENLARGE - enlarge the axes of a figure (f,h)
- v_randfilt V_RANDFILT Generate filtered gaussian noise without initial transient
- v_spgrambw V_SPGRAMBW Draw spectrogram [T,F,B]=(s,fs,mode,bw,fmax,db,tinc,ann)
- v_texthvc V_TEXTHVC - write text on graph with specified alignment and colour
- v_xticksi V_XTIXKSI labels the x-axis of a plot using SI multipliers S=(AH)

This function is called by:

- v_activlev V_ACTIVLEV Measure active speech level as in ITU-T P.56 [LEV,AF,FSO]=(sp,FS,MODE)
- v_activlevg V_ACTIVLEVG Measure active speech level robustly [LEV,AF,FSO]=(sp,FS,MODE)
- v_addnoise V_ADDNOISE Add noise at a chosen SNR [z,p,fso]=(s,fsx,snr,m,nb,fsa)
- v_earnoise V_EARNOISE Add noise to simulate the hearing threshold of a listener [Y,X,V]=(S,FS,M,SPL)
- v_fxpefac V_FXPEFAC PEFAC pitch tracker [FX,TT,PV,FV]=(S,FS,TINC,M,PP)
- v_sigalign V_SIGALIGN align a clean reference with a noisy signal [d,g,rr,ss]=(s,r,maxd,m,fs)
- v_snrseg V_SNRSEG Measure segmental and global SNR [SEG,GLO]=(S,R,FS,M,TF)

## 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 Lesser General Public License as published by
    0124 %   the Free Software Foundation; either version 3 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 Lesser General Public License for more details.
    0131 %
    0132 %   You can obtain a copy of the GNU Lesser General Public License from
    0133 %   https://www.gnu.org/licenses/ .
    0134 %
    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