v_randfilt

PURPOSE ^

V_RANDFILT Generate filtered gaussian noise without initial transient

SYNOPSIS ^

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

DESCRIPTION ^

V_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() as given below.
 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.

 The state space representation of filter() is z=Az+Bx, y=Cz+Dx where
       A=[-pa(2:k); eye(k-2) zeros(k-2,1)].'
       B=[pb(2:k)-pb(1)*pa(2:k)].'
       C=[1 zeros(1,k-2)]
       D=pb(1)
 and k is the order of the filter, pa(1)=1 and pa and pb are zero-padded to length k+1.

 Refs:
 [1]    D. M. Brookes. Coloured noise generation without filter startup transient.
       IEE Electronics Lett., 37 (4): 255256, Feb. 2001. doi: 10.1049/el:20010144.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [y,zf,u,p]=v_randfilt(pb,pa,ny,zi)
0002 %V_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() as given below.
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 % The state space representation of filter() is z=Az+Bx, y=Cz+Dx where
0020 %       A=[-pa(2:k); eye(k-2) zeros(k-2,1)].'
0021 %       B=[pb(2:k)-pb(1)*pa(2:k)].'
0022 %       C=[1 zeros(1,k-2)]
0023 %       D=pb(1)
0024 % and k is the order of the filter, pa(1)=1 and pa and pb are zero-padded to length k+1.
0025 %
0026 % Refs:
0027 % [1]    D. M. Brookes. Coloured noise generation without filter startup transient.
0028 %       IEE Electronics Lett., 37 (4): 255256, Feb. 2001. doi: 10.1049/el:20010144.
0029 
0030 %      Copyright (C) Mike Brookes 1997-2021
0031 %      Version: $Id: v_randfilt.m 10865 2018-09-21 17:22:45Z dmb $
0032 %
0033 %   VOICEBOX is a MATLAB toolbox for speech processing.
0034 %   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
0035 %
0036 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0037 %   This program is free software; you can redistribute it and/or modify
0038 %   it under the terms of the GNU General Public License as published by
0039 %   the Free Software Foundation; either version 2 of the License, or
0040 %   (at your option) any later version.
0041 %
0042 %   This program is distributed in the hope that it will be useful,
0043 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0044 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0045 %   GNU General Public License for more details.
0046 %
0047 %   You can obtain a copy of the GNU General Public License from
0048 %   http://www.gnu.org/copyleft/gpl.html or by writing to
0049 %   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
0050 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0051 
0052 % first normalize the denominator polynomial if necessary
0053 
0054 if pa(1)~=1
0055     pb=pb/pa(1); pa=pa/pa(1);
0056 end
0057 
0058 % check to see if we must generate zi
0059 
0060 if nargin<4 | nargout>2
0061     lb=length(pb);
0062     la=length(pa);
0063 
0064     k=max(la,lb)-1;     % filter order
0065     n=la-1;             % denominator order
0066     ii=k+1-n:k;
0067 
0068     % form controllability matrix
0069 
0070     q=zeros(k,k);
0071     [z,q(:,1)]=filter(pb,pa,1);
0072     for i=2:k [z,q(:,i)]=filter(pb,pa,0,q(:,i-1)); end
0073 
0074     % we generate m through the step-down procedure
0075     s=randn(k,1);
0076     if n
0077         m=zeros(n,n);
0078         g=pa;
0079         for i=1:n
0080             g=(g(1)*g(1:end-1)-g(end)*g(end:-1:2))/sqrt((g(1)-g(end))*(g(1)+g(end)));
0081             m(i,i:n)=g;
0082         end
0083         s(ii)=triu(toeplitz(pa(1:n)))*(m\s(ii));
0084         if nargout>2
0085             u=q;
0086             u(:,ii)=q(:,ii)*triu(toeplitz(pa(1:n)))/m;
0087         end
0088     else
0089         if nargout>2
0090             u=q;
0091         end
0092     end
0093     if nargin < 4
0094         if k
0095             zi=q*s;
0096         else
0097             zi=[];
0098         end
0099     end
0100 end
0101 if nargout>2
0102     if ~numel(u)
0103         p=pb(1).^2;
0104     else
0105         p=u(1,:)*u(1,:)'+pb(1).^2;
0106     end
0107 end
0108 if nargin>2 && ny>0
0109     [y,zf]=filter(pb,pa,randn(ny,1),zi);
0110 else
0111     zf=zi;
0112     y=[];
0113 end

Generated by m2html © 2003