V_RANDVEC Generate real or complex GMM/lognormal random vectors X=(N,M,C,W,MODE) generates a random matrix of size (|n|,p) where p is the maximum dimension of M or C (see note below about row versus column vectors) Inputs: N is the number of points to generate M(K,P) is the mean vectors (one row per mixture) C(K,P) are diagonal covariances (one row per mixture) or C(P,P,K) are full covariance matrices (one per mixture) W(K) are the mixture weights (or omit if all mixtures have equal weight) MODE character string specifying options: g = real gaussian [default] c = complex gaussian l = lognormal Outputs: X(N,P) is the output data KX(N,1) Mixture number associated with each row of X Note this routine generates row vectors such that E((x-m)'*(x-m)) = C = cov(x). If Alternatively x=(n,m',c,w,mode)' will generate column vectors satisfying E((x-m)*(x-m)') = C = cov(x').
0001 function [x,kx]=v_randvec(n,m,c,w,mode) 0002 %V_RANDVEC Generate real or complex GMM/lognormal random vectors X=(N,M,C,W,MODE) 0003 % generates a random matrix of size (|n|,p) where p is the maximum 0004 % dimension of M or C (see note below about row versus column vectors) 0005 % Inputs: N is the number of points to generate 0006 % M(K,P) is the mean vectors (one row per mixture) 0007 % C(K,P) are diagonal covariances (one row per mixture) 0008 % or C(P,P,K) are full covariance matrices (one per mixture) 0009 % W(K) are the mixture weights (or omit if all mixtures have equal weight) 0010 % MODE character string specifying options: 0011 % g = real gaussian [default] 0012 % c = complex gaussian 0013 % l = lognormal 0014 % 0015 % Outputs: X(N,P) is the output data 0016 % KX(N,1) Mixture number associated with each row of X 0017 % 0018 % Note this routine generates row vectors such that E((x-m)'*(x-m)) = C = cov(x). If 0019 % Alternatively x=(n,m',c,w,mode)' will generate column vectors satisfying E((x-m)*(x-m)') = C = cov(x'). 0020 0021 % Bugs/suggestions 0022 % (1) New mode 'x' to approximate chi squared 0023 0024 % Copyright (C) Mike Brookes 1998 0025 % Version: $Id: v_randvec.m 10865 2018-09-21 17:22:45Z dmb $ 0026 % 0027 % VOICEBOX is a MATLAB toolbox for speech processing. 0028 % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 0029 % 0030 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0031 % This program is free software; you can redistribute it and/or modify 0032 % it under the terms of the GNU General Public License as published by 0033 % the Free Software Foundation; either version 2 of the License, or 0034 % (at your option) any later version. 0035 % 0036 % This program is distributed in the hope that it will be useful, 0037 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0038 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0039 % GNU General Public License for more details. 0040 % 0041 % You can obtain a copy of the GNU General Public License from 0042 % http://www.gnu.org/copyleft/gpl.html or by writing to 0043 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0044 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0045 0046 % first sort out the input arguments 0047 0048 sm=size(m); 0049 if nargin<3 0050 c=ones(sm); % default to unit variance 0051 end 0052 sc=size(c); 0053 p=max(sm(2),sc(2)); % data dimension 0054 k=sm(1); % number of mixtures 0055 fullc=(length(sc)>2) || (sc(1)>k); 0056 if nargin<4 0057 mode='g'; % default to gaussian 0058 w=ones(k,1); 0059 else 0060 if ischar(w) 0061 mode = w; % w argument has been omitted 0062 w=ones(k,1); 0063 elseif nargin<5 0064 mode='g'; 0065 end 0066 end 0067 ty=mode(1); % ignore all but first character for type 0068 x=zeros(n,p); % initialize output array 0069 if sm(2)~=p 0070 m=repmat(m,1,p); % if m is given as a scalar 0071 end 0072 if sc(2) ~=p 0073 c=repmat(c,1,p); % if c is given as a scalar 0074 end 0075 q=sqrt(0.5); 0076 if k>1 0077 kx=v_randiscr(w,n); 0078 else 0079 kx=ones(n,1); 0080 end 0081 for kk=1:k 0082 nx=find(kx==kk); 0083 nn=length(nx); 0084 if nn % check if we need to generate any from mixture kk 0085 0086 % extract the mean and cov for this mixture 0087 0088 mm=m(kk,:); % mean vector 0089 if fullc % full covariance matrix 0090 cc=c(:,:,kk); 0091 if ty=='l' % lognormal distribution - convert mean and covariance 0092 cc=log(1+cc./(mm.'*mm)); 0093 mm=log(mm)-0.5*diag(cc).'; 0094 end 0095 else 0096 cc=c(kk,:); 0097 if ty=='l' % lognormal distribution - convert mean and covariance 0098 cc=log(1+cc(:).'./mm.^2); 0099 mm=log(mm)-0.5*cc; 0100 end 0101 end 0102 0103 % now generate nn complex or real values 0104 0105 if ty=='c' % generate complex values 0106 xx=q*randn(nn,p)+1i*q*randn(nn,p); % complex-valued unit variance values 0107 else 0108 xx=randn(nn,p); % real-valued unit variance values 0109 end; 0110 0111 % scale by the square root of the covariance matrix 0112 0113 if fullc % full covariance covariance 0114 [v,d]=eig((cc+cc')/2); % force covariance matrix to be hermitian 0115 xx=(xx.*repmat(sqrt(abs(diag(d))).',nn,1))*v'+repmat(mm,nn,1); % and also positive definite 0116 else 0117 xx=xx.*repmat(sqrt(abs(cc)),nn,1)+repmat(mm,nn,1); % different mean/cov for each column 0118 end 0119 x(nx,:)=xx; 0120 end 0121 end 0122 if ty=='l' % lognormal distribution 0123 x=exp(x); 0124 end 0125 if ~nargout 0126 if p==1 0127 if ty=='c' 0128 plot(real(x), imag(x),'+'); 0129 xlabel('Real'); 0130 ylabel('Imag'); 0131 else 0132 nbin=max(min(floor(sqrt(n)),50),5); 0133 hist(x,nbin); 0134 xlabel('X'); 0135 ylabel('Frequency'); 0136 end 0137 else 0138 [vv,iv]=sort(var(x,0,1)); 0139 iv=sort(iv(end-1:end)); 0140 plot(real(x(:,iv(1))), real(x(:,iv(2))),'+'); 0141 if ty=='c' 0142 xylab='Real[ x(%d) ]'; 0143 else 0144 xylab='x(%d)'; 0145 end 0146 xlabel(sprintf(xylab,iv(1))); 0147 ylabel(sprintf(xylab,iv(2))); 0148 end 0149 end 0150