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
                    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
                    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
                    A - plot zeros/poles
                    U - plot using a uniform frequency axis
                    L - plot using a log frequency axis
                    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].
                     Converted from mouth reference point @ 0.025m to dB SPL @ 1m on-axis.
  13  LTASS-1994   : the long-term average speech spectrum that is taken from Table 2 in [2]

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

Generated on Tue 10-Oct-2017 08:30:10 by m2html © 2003