# 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): 255–256, Feb. 2001. doi: 10.1049/el:20010144.

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

- v_addnoise V_ADDNOISE Add noise at a chosen SNR [z,p,fso]=(s,fsx,snr,m,nb,fsa)
- v_stdspectrum V_STDSPECTRUM Generate standard acoustic/speech spectra in s- or z-domain [B,A,SI,SN]=(S,M,F,N,ZI,BS,AS)
- v_usasi V_USASI generates N samples of USASI noise at sample frequency FS X=(N,FS)

## 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): 255–256, 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 Lesser General Public License as published by
    0039 %   the Free Software Foundation; either version 3 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 Lesser General Public License for more details.
    0046 %
    0047 %   You can obtain a copy of the GNU Lesser General Public License from
    0048 %   https://www.gnu.org/licenses/ .
    0049 %
    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