[H,D,F] = GenerateMIMOSources(L,P,M,K,gamma,Mode)
GenerateMIMOSources(L,P,M,K) produces a source model for M-array data
emanating from L independent Gaussian sources spectrally shaped by Pth
order moving average (MA, i.e. finite impulse response) innovation
filters. The source signals are then mixed by a random MxL paraunitary
matrix of order K.
GenerateMIMOSources(L,P,M,K,gamma) operates the same but the additional
input parameter gamma limits the radii of zeros in the MA innovation
filters. This restricts the dynamic range of the generated spectra.
A flow graph of this source model is given in [1]. In the context of
multichannel coding, M=L is required to avoid a division by zero, as the
geometric mean of variances is zero for M>L. The purpose of limiting the
dynamic range of source PSDs is driven by the same aim, and keeps the
problem numerically tractable.
The source model does currently not permit the overdetermined case where
the number of sources exceeds the number of sensors.
[H,D,F]=GenerateMIMOSources(L,P,M,K,gamma) returns the paraunitary
mixing matrix H(z) of dimension MxL and order K in H; D represents a
diagonal parahermitian matrix D(z) containing the power spectra of the L
sources prior to convolutive mixing. Minimum phase filters of order P
that generate such PSDs from unit variance complex Gaussian white noise
are contained in the columns of the matrix F.
[H,D,F]=GenerateMIMOSources(L,P,M,K,gamma,Mode) returns a real-valued
source model for Mode='real', and a complex valued one by default.
Input parameters
L number of independent sources
P order of MA innovation filter
M dimension of array generated (M>=L)
K order of paraunitary mixing matrix
gamma max radius for zeros in the innovation filters
(default value is 1)
Mode 'real' for readl valued, 'complex valued' (default)
Output parameters
H paraunitary matrix
D diagonal, parahermitian matrix
F columns are the innovation filters, whose
ACFs generate the diagonal entries of D.
Reference:
[1] S. Redif, S. Weiss and J.G. McWhirter, "Sequential Matrix
Diagonalisation Algorithms for Polynomial EVD of Parahermitian
Matrices," IEEE Transactions on Signal Processing, 63(1):81-89,
January 2015.
0001 function [H,D,F] = GenerateMIMOSources(L,P,M,K,gamma,Mode) 0002 %[H,D,F] = GenerateMIMOSources(L,P,M,K,gamma,Mode) 0003 % 0004 % GenerateMIMOSources(L,P,M,K) produces a source model for M-array data 0005 % emanating from L independent Gaussian sources spectrally shaped by Pth 0006 % order moving average (MA, i.e. finite impulse response) innovation 0007 % filters. The source signals are then mixed by a random MxL paraunitary 0008 % matrix of order K. 0009 % 0010 % GenerateMIMOSources(L,P,M,K,gamma) operates the same but the additional 0011 % input parameter gamma limits the radii of zeros in the MA innovation 0012 % filters. This restricts the dynamic range of the generated spectra. 0013 % 0014 % A flow graph of this source model is given in [1]. In the context of 0015 % multichannel coding, M=L is required to avoid a division by zero, as the 0016 % geometric mean of variances is zero for M>L. The purpose of limiting the 0017 % dynamic range of source PSDs is driven by the same aim, and keeps the 0018 % problem numerically tractable. 0019 % 0020 % The source model does currently not permit the overdetermined case where 0021 % the number of sources exceeds the number of sensors. 0022 % 0023 % [H,D,F]=GenerateMIMOSources(L,P,M,K,gamma) returns the paraunitary 0024 % mixing matrix H(z) of dimension MxL and order K in H; D represents a 0025 % diagonal parahermitian matrix D(z) containing the power spectra of the L 0026 % sources prior to convolutive mixing. Minimum phase filters of order P 0027 % that generate such PSDs from unit variance complex Gaussian white noise 0028 % are contained in the columns of the matrix F. 0029 % 0030 % [H,D,F]=GenerateMIMOSources(L,P,M,K,gamma,Mode) returns a real-valued 0031 % source model for Mode='real', and a complex valued one by default. 0032 % 0033 % Input parameters 0034 % L number of independent sources 0035 % P order of MA innovation filter 0036 % M dimension of array generated (M>=L) 0037 % K order of paraunitary mixing matrix 0038 % gamma max radius for zeros in the innovation filters 0039 % (default value is 1) 0040 % Mode 'real' for readl valued, 'complex valued' (default) 0041 % 0042 % Output parameters 0043 % H paraunitary matrix 0044 % D diagonal, parahermitian matrix 0045 % F columns are the innovation filters, whose 0046 % ACFs generate the diagonal entries of D. 0047 % 0048 % Reference: 0049 % 0050 % [1] S. Redif, S. Weiss and J.G. McWhirter, "Sequential Matrix 0051 % Diagonalisation Algorithms for Polynomial EVD of Parahermitian 0052 % Matrices," IEEE Transactions on Signal Processing, 63(1):81-89, 0053 % January 2015. 0054 0055 % Stephan Weiss, August 2012 0056 % real valued option, S. Weiss, Oct 2018 0057 0058 %-------------------------------------------------------- 0059 % Input Parameter Check 0060 %-------------------------------------------------------- 0061 if nargin==4, 0062 gamma=1; 0063 end; 0064 if M < L, 0065 error('more sources than array elements'); 0066 end; 0067 RealMode = 0; 0068 if nargin==6, 0069 if Mode=='real', 0070 RealMode=1; 0071 end; 0072 end; 0073 0074 %-------------------------------------------------------- 0075 % Paraunitary Matrix 0076 %-------------------------------------------------------- 0077 H = eye(M); 0078 for k = 1:K, 0079 if RealMode==1, 0080 v = randn(M,1); 0081 else 0082 v = randn(M,1) + sqrt(-1)*randn(M,1); 0083 end; 0084 v = v./norm(v,2); 0085 Hnew = zeros(M,M,2); 0086 if RealMode==1, 0087 [Q,~] = qr(randn(M,M)); 0088 else 0089 [Q,~] = qr(randn(M,M) + sqrt(-1)*randn(M,M)); 0090 end; 0091 Hnew(:,:,1) = v*(v'*Q); 0092 Hnew(:,:,2) = Q - Hnew(:,:,1); 0093 H = PolyMatConv(H,Hnew); 0094 end; 0095 0096 %-------------------------------------------------------- 0097 % Source Innovation Filters 0098 %-------------------------------------------------------- 0099 if RealMode==1, 0100 P2 = floor(P/2); 0101 P = 2*P2; 0102 Fz = (rand(P2-1,L)*gamma).*exp(sqrt(-1)*2*pi*rand(P2-1,L)); % zeros 0103 Fz = [Fz; conj(Fz)]; 0104 else 0105 Fz = (rand(P-1,L)*gamma).*exp(sqrt(-1)*2*pi*rand(P-1,L)); % zeros 0106 end; 0107 F = zeros(P,L); 0108 for l = 1:L, 0109 F(:,l) = poly( Fz(:,l)).'; 0110 end; 0111 F(:,1) = F(:,1)./norm(F(:,1),2); 0112 Ffd = abs(fft(F,1024,1)); 0113 for l = 1:L-1, 0114 alpha = min(Ffd(:,l)./Ffd(:,l+1)); 0115 Ffd(:,l+1) = alpha*Ffd(:,l+1); 0116 F(:,l+1) = alpha*F(:,l+1); 0117 end; 0118 0119 %-------------------------------------------------------- 0120 % Output Parameters 0121 %-------------------------------------------------------- 0122 D1 = zeros(L,L,P); 0123 for l = 1:L, 0124 for p = 1:P, 0125 D1(l,l,p) = F(p,l); 0126 end; 0127 end; 0128 D = PolyMatConv(D1,ParaHerm(D1)); 0129 H = H(:,1:L,:); 0130 0131 0132