Reply by Rune Allnor June 17, 20092009-06-17
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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%