Home > voicebox > randfilt.m

# randfilt

## PURPOSE

RANDFILT Generate filtered gaussian noise without initial transient

## SYNOPSIS

function [y,zf,u,p]=randfilt(pb,pa,ny,zi)

## DESCRIPTION

```RANDFILT Generate filtered gaussian noise without initial transient

Inputs: pb(1,:)  Numerator polynomial of discrete time filter
pa(1,:)  Denominator polynomial of discrete time filter
ny       Number of output samples required
zi       Optional initial filter state

Outputs: y(ny,1)  Filtered Gaussian noise
zf       final filter state (can be used to extend the noise sequence)
u        The state covariance matrix, <zf*zf'>, is u*u'
p        Is the expected value of y(i)^2

zf and zi are the output and optional input state as defined in filter()
If zi is not specified, random numbers with the correct covariance will be used.
output u*u' is the state covariance matrix for filter(). Output p is the
mean power of the output signal y.```

## CROSS-REFERENCE INFORMATION

This function calls:
This function is called by:
• stdspectrum STDSPECTRUM Generate standard acoustic/speech spectra in s- or z-domain [B,A,SI,SN]=(S,M,F,N,ZI,BS,AS)
• usasi USASI generates N samples of USASI noise at sample frequency FS X=(N,FS)

## SOURCE CODE

```0001 function [y,zf,u,p]=randfilt(pb,pa,ny,zi)
0002 %RANDFILT Generate filtered gaussian noise without initial transient
0003 %
0004 %  Inputs: pb(1,:)  Numerator polynomial of discrete time filter
0005 %          pa(1,:)  Denominator polynomial of discrete time filter
0006 %          ny       Number of output samples required
0007 %          zi       Optional initial filter state
0008 %
0009 % Outputs: y(ny,1)  Filtered Gaussian noise
0010 %          zf       final filter state (can be used to extend the noise sequence)
0011 %          u        The state covariance matrix, <zf*zf'>, is u*u'
0012 %          p        Is the expected value of y(i)^2
0013 %
0014 % zf and zi are the output and optional input state as defined in filter()
0015 % If zi is not specified, random numbers with the correct covariance will be used.
0016 % output u*u' is the state covariance matrix for filter(). Output p is the
0017 % mean power of the output signal y.
0018
0019 %      Copyright (C) Mike Brookes 1997
0020 %      Version: \$Id: randfilt.m 5885 2015-03-18 09:35:01Z dmb \$
0021 %
0022 %   VOICEBOX is a MATLAB toolbox for speech processing.
0024 %
0025 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0026 %   This program is free software; you can redistribute it and/or modify
0028 %   the Free Software Foundation; either version 2 of the License, or
0029 %   (at your option) any later version.
0030 %
0031 %   This program is distributed in the hope that it will be useful,
0032 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0033 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0034 %   GNU General Public License for more details.
0035 %
0036 %   You can obtain a copy of the GNU General Public License from
0037 %   http://www.gnu.org/copyleft/gpl.html or by writing to
0038 %   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
0039 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0040
0041 % first normalize the denominator polynomial if necessary
0042
0043 if pa(1)~=1
0044     pb=pb/pa(1); pa=pa/pa(1);
0045 end
0046
0047 % check to see if we must generate zi
0048
0049 if nargin<4 | nargout>2
0050     lb=length(pb);
0051     la=length(pa);
0052
0053     k=max(la,lb)-1;
0054     l=la-1;
0055     ii=k+1-l:k;
0056
0057     % form controllability matrix
0058
0059     q=zeros(k,k);
0060     [z,q(:,1)]=filter(pb,pa,1);
0061     for i=2:k [z,q(:,i)]=filter(pb,pa,0,q(:,i-1)); end
0062
0063     % we generate m through the step-down procedure
0064     s=randn(k,1);
0065     if l
0066         m=zeros(l,l);
0067         g=pa;
0068         for i=1:l
0069             g=(g(1)*g(1:end-1)-g(end)*g(end:-1:2))/sqrt((g(1)-g(end))*(g(1)+g(end)));
0070             m(i,i:l)=g;
0071         end
0072         s(ii)=triu(toeplitz(pa(1:l)))*(m\s(ii));
0073         if nargout>2
0074             u=q;
0075             u(:,ii)=q(:,ii)*triu(toeplitz(pa(1:l)))/m;
0076         end
0077     else
0078         if nargout>2
0079             u=q;
0080         end
0081     end
0082     if nargin < 4
0083         if k
0084             zi=q*s;
0085         else
0086             zi=[];
0087         end
0088     end
0089 end
0090 if nargout>2
0091     if ~numel(u)
0092         p=pb(1).^2;
0093     else
0094         p=u(1,:)*u(1,:)'+pb(1).^2;
0095     end
0096 end
0097 if nargin>2 && ny>0
0098     [y,zf]=filter(pb,pa,randn(ny,1),zi);
0099 else
0100     zf=zi;
0101     y=[];
0102 end```

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