p_music.m

现代数字信号处理中基于互相关函数的music方法

clc
clear
var = sqrt(1/exp(1.0));
v = var*randn(1,1000);

n = 1:1000;
w0 = 1.4*pi;
xn = exp(j*pi*n-j*pi)+exp(j*w0*n-j*0.7*pi)+2*exp(j*1.8*n-j*0.7*pi)+v;
% mm1=[0.697 0.68 0.71 0.694];
% xn=10*exp(-0.0027*n).*cos(2*pi*mm1(1)*n+pi/5)+2*exp(-0.2652*n).*cos(2*pi*mm1(2)*n+pi/6)+10*exp(-0.0311*n).*cos(2*pi*mm1(3)*n+pi/7)+exp(-0.2936*n).*cos(2*pi*mm1(4)*n+pi/4)+v;
m = -500:500;
[X,R] = corrmtx(xn,20,'covariance');
% [r,lag] = xcorr(xn,500,'biased');
% 
% %pmusic
% R = [r(501),r(502),r(503),r(504),;
%      r(500),r(501),r(502),r(503);
%      r(499),r(500),r(501),r(502);
%      r(498),r(499),r(500),r(501)];
 
[U,V] = eig(R);
% V3 = V(:,3);
% V4 = V(:,4);
Ug = U(:,4:end);

p = 0:20;
wm = [0:0.002*pi:2*pi];
B = [(exp(-j)).^(wm'*p)];
A = B.';
toc;

ss=0;
for k = 1:18
    p1 = A'*U(:,k+1);
    ss = ss + abs(p1).^2;
    k = k+1;
end
pmus = 1./ss';

subplot(211);
pp = 10*log10(pmus);
plot(wm/pi);
ylabel('Power Spectrum Magnitude(dB)');
title('Music Spectrum Estimate');

%
subplot(212);
pmusic(xn,20);