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.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 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