PssMusic

 function [Pss,Ps] = PssMusic(Q,L,K,angles);

 Polynomial MUSIC algorithm.

 Input parameter:
   Q      polynomial (spatio-temporal) modal matrix
   L      number of sources
   K      number of frequency bins evaluated
   angles possible angles of arrival to be evaluated
          (measured in rad against broadside)
   epsilon   regularisation constant to avoid division by zero

 Output parameters:
   Pss    polynomial spatio-spectral music
   Ps     polynomial spatial music

 S. Weiss, University of Strathclyde, 30/1/2013

0001 function [Pss,Ps,AA] = PssMusic(Q,L,K,angles,epsilon);
0002 % function [Pss,Ps] = PssMusic(Q,L,K,angles);
0003 %
0004 % Polynomial MUSIC algorithm.
0005 %
0006 % Input parameter:
0007 %   Q      polynomial (spatio-temporal) modal matrix
0008 %   L      number of sources
0009 %   K      number of frequency bins evaluated
0010 %   angles possible angles of arrival to be evaluated
0011 %          (measured in rad against broadside)
0012 %   epsilon   regularisation constant to avoid division by zero
0013 %
0014 % Output parameters:
0015 %   Pss    polynomial spatio-spectral music
0016 %   Ps     polynomial spatial music
0017 %
0018 % S. Weiss, University of Strathclyde, 30/1/2013
0019 
0020 M = size(Q,1);
0021 
0022 Qn = Q(L+1:M,:,:);
0023 Rn = MIMOConv(ParaHerm(Qn),Qn);
0024 
0025 
0026 for i = 1:size(angles,2),
0027     X = BBSteerVec(M,angles(i),100);
0028     Y =  MIMOConv(ParaHerm(X),MIMOConv(Rn,X));
0029     dummy2 = shiftdim(Y,1);
0030     AA(i,:) = dummy2;
0031 end;
0032 
0033 if K < size(AA,2),
0034    disp(['warning: FFT length is shorter than auto-correlation sequence ' ...
0035          'in function PssMusic()']);
0036 end;
0037 
0038 Pss = 1./(abs(fft(AA,K,2).') + epsilon);
0039 Ps = 1./abs(AA(:,(size(AA,2)+1)/2) + epsilon);
0040 
0041 
0042 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0043 function X = BBSteerVec(M,alpha,T);
0044 X = zeros(M,1,2*T+1);
0045 dT= sin(alpha)/2;
0046 t = (-T:T);
0047 Hann = .5 + .5*cos(pi*t/T);
0048 TotalNorm = 0;
0049 for m = 1:M,
0050   dummy = SincFct(t-(m-1)*dT);  %.*Hann;
0051   X(m,1,:) = dummy;
0052   %  TotalNorm = TotalNorm + norm(dummy,2);
0053 end;
0054 %X = X./TotalNorm;