# How MUSIC & friends fail

Started by June 17, 2009
Hi all.

Over the past few days I have commented on MUSIC and its
pathologies in at least two dofferent threads.  I have suggested
that people do a certain excercise to see this pathology for
themselves.

While this is useful, it requires that the users

1) Know what to look for
2) Actually do the job of setting up the experiment

So with that in mind, I have set up a very simple matlab demo.

The demo does only require a basic lisence (no SP toolboxes),
and I have made the effort to keep it as basic as possible, so that
the thing runs on as many matlab versions as possible. Save the
functions to separate files, or use one file with nested functions.

Just make sure that the main function is MUSICFailDemo.

That said, I'm using R2006a, so it might be that users of earlier
versions might get problems.

The demo is designed as follows:

1) The main function is MUSICFailDemo, where the user chooses
the simulation parameters

P - order of covariance matrix
N - Number of samples in data series
sigma2 - Variance of noise.

One can choose to run the simulation with default parameters.

2) The main loop adds succesively more signal components
to the signal, and runs the simulation with 50 different
realizations of noise. (Even if sigma2 = 0).

3) The resulting frequency estimates are plotted, along with a grid
where the expected frequency estimates:

- The horizontal green lines show where frequency estimates
should end up in a perfect world
- The black vertical lines show the transitions between numbers
of signal components
- The red vertical line shows where the complex-valued signal
violates the assumptions/conditions in the MUSIC model
- The blue dots are the actual frequency estimates.

Note also that I have skipped the order estimation problem:
In the simulation, the MUSIC estimator is set to look for the
true number of signal components.

It's a rather basic demo that ought to be

1) Repeated with 'proper' implementations of MUSIC, ESPRIT etc.
2) Included in the textbooks.

Have fun.

Rune

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function MUSICFailDemo(P,N,sigma2)
% MUSICFailDemo - Demonstrates why MUSIC is
% dangerous in real-world applications.
% P sinusoidals are successively added to the
% signal. Each sinusoidal has unit amplitude.
%
% P      - Order of covariance matrix (Default = 8)
% N      - Number of samples in data sequence (Default = 4*P)
% sigma2 - Variance of Gaussian noise (Default = 0)

if nargin < 1
P = 8;
else
if ~(isnumeric(P) & sum((size(P) == [1,1]))==2 & P > 0)
error('P must be positive scalar')
end
end

if nargin < 2
N = 4*P;
else
if ~(isnumeric(N) & sum((size(N) == [1,1]))==2 & N > 0)
error('N must be positive scalar')
end
end

if nargin < 3
sigma2 = 0;
else
if ~(isnumeric(sigma2) & sum((size(sigma2) == [1,1]))==2 &
sigma2 >= 0)
error('sigma2 must be non-negative scalar')
end
end

Q = 50;
Wm = zeros(2*P,Q*2*P);
wvec = 2*pi*(1:(2*P))/(2*P+1);
r = 1;
for p = 1:(2*P)
for q = 1:Q
s = makesignal(wvec(1:p),N)+sigma2*randn(N,1);
wv=publicMUSIC(s,p,P,8192);
Wm(1:min(p,length(wv)),(p-1)*Q+q) = wv;
end
end
cla
for p=1:(2*P)
pidx = find(Wm(p,:)~=0);
plot(pidx/Q,Wm(p,pidx),'.')
hold on
plot([(p-1),2*P],wvec(p)*[1,1],'g','linewidth',2)
end
ax = axis;
for p=0:(2*P)
plot(p*[1,1],ax(3:4),'k')
end
plot((P-1)*[1,1],ax(3:4),'--r','linewidth',3)
title(['P = ',int2str(P),' sigma2 = ',num2str(sigma2)])
xlabel('# of signal components')
ylabel('Frequency $[0,2\pi]$','interpreter','latex')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function s = makesignal(wv,N)
s = exp(j*reshape(0:(N-1),N,1)*reshape(wv,1,length(wv)))*ones
(length(wv),1);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function rxx = covarmat(x,P)
rxx = zeros(length(x)+P-1,P);
for n=1:(P)
rxx(:,n)=[zeros(n-1,1);x;zeros(P-n,1)];
end
rxx=rxx(P:end-P+1,:)'*rxx(P:end-P+1,:);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function wv = findminima(Nv,N)
[Mn,Nn]=size(Nv);
Iv = zeros(N+2,1);
for n=1:N
Iv(n+1) = abs(sum(Nv*exp(j*2*pi*(reshape(0:(Nn-1),Nn,1)*(n-1)/
N))));
end
Iv(1)=Iv(end-1);
Iv(end) = Iv(2);
midx = find(Iv(2:end-1)<Iv(1:end-2) & Iv(2:end-1) < Iv(3:end));
Mv = Iv(midx);
[mm,smidx]=sort(Mv);

wv = 2*pi*(midx(smidx)-1)/N;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function wv = publicMUSIC(x,d,P,Nwv)
x = reshape(x,length(x),1);
rxx = covarmat(x,P);
[U,S,V] = svd(rxx);
Nv = (V(:,d+1:end))';
wv = findminima(Nv,Nwv);
wv=wv(1:min(d,length(wv)));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%