RANDFILT Generate filtered gaussian noise without initial transient Generates noise by passing gaussian noise through a filter pb/pa output signal y is a column vector of length ny 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().
0001 function [y,zf,u]=randfilt(pb,pa,ny,zi) 0002 %RANDFILT Generate filtered gaussian noise without initial transient 0003 % 0004 % Generates noise by passing gaussian noise through a filter pb/pa 0005 % output signal y is a column vector of length ny 0006 % 0007 % zf and zi are the output and optional input state as defined in filter() 0008 % If zi is not specified, random numbers with the correct covariance will be used. 0009 % output u*u' is the state covariance matrix for filter(). 0010 0011 % Copyright (C) Mike Brookes 1997 0012 % Version: $Id: randfilt.m 713 2011-10-16 14:45:43Z dmb $ 0013 % 0014 % VOICEBOX is a MATLAB toolbox for speech processing. 0015 % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 0016 % 0017 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0018 % This program is free software; you can redistribute it and/or modify 0019 % it under the terms of the GNU General Public License as published by 0020 % the Free Software Foundation; either version 2 of the License, or 0021 % (at your option) any later version. 0022 % 0023 % This program is distributed in the hope that it will be useful, 0024 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0025 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0026 % GNU General Public License for more details. 0027 % 0028 % You can obtain a copy of the GNU General Public License from 0029 % http://www.gnu.org/copyleft/gpl.html or by writing to 0030 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0031 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0032 0033 % check to see if we must generate zi 0034 0035 if pa(1)~=1 0036 pb=pb/pa(1); pa=pa/pa(1); 0037 end 0038 if nargin<4 | nargout>2 0039 lb=length(pb); 0040 la=length(pa); 0041 0042 k=max(la,lb)-1; 0043 l=la-1; 0044 ii=k+1-l:k; 0045 0046 % form controllability matrix 0047 0048 q=zeros(k,k); 0049 [z,q(:,1)]=filter(pb,pa,1); 0050 for i=2:k [z,q(:,i)]=filter(pb,pa,0,q(:,i-1)); end 0051 0052 % we generate m through the step-down procedure 0053 s=randn(k,1); 0054 if l 0055 m=zeros(l,l); 0056 g=pa; 0057 for i=1:l 0058 g=(g(1)*g(1:end-1)-g(end)*g(end:-1:2))/sqrt((g(1)-g(end))*(g(1)+g(end))); 0059 m(i,i:l)=g; 0060 end 0061 s(ii)=triu(toeplitz(pa(1:l)))*(m\s(ii)); 0062 if nargout>2 0063 u=q; 0064 u(:,ii)=q(:,ii)*triu(toeplitz(pa(1:l)))/m; 0065 end 0066 else 0067 if nargout>2 0068 u=q; 0069 end 0070 end 0071 if nargin < 4 0072 if k 0073 zi=q*s; 0074 else 0075 zi=[]; 0076 end 0077 end 0078 end 0079 if ny 0080 [y,zf]=filter(pb,pa,randn(ny,1),zi); 0081 else 0082 zf=zi; 0083 y=[]; 0084 end