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:

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.
0023 %   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
0024 %
0025 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0026 %   This program is free software; you can redistribute it and/or modify
0027 %   it under the terms of the GNU General Public License as published by
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 Tue 10-Oct-2017 08:30:10 by m2html © 2003