Home > voicebox > stdspectrum.m

stdspectrum

PURPOSE ^

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]=stdspectrum(s,m,f,n,zi,bs,as)

DESCRIPTION ^

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

Generated on Mon 06-Aug-2018 14:48:32 by m2html © 2003