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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%