music.m

本文提供一种用MUSIC算法恢复谐波的一个实例

clear all
n=1:128;
N1=128;
m=30;
%p=2;
N=50;
%设置观测数据
x=sqrt(20)*sin(2*pi*0.2*n)+sqrt(2)*sin(2*pi*0.213*n)+randn(1,N1);
%建立自相关矩阵
Rxx=zeros(m);
for k=1:N;
        for i=1:m;
             xm(i)=x(i+k-1);
        end
    rxx=xm'*xm;
    Rxx=Rxx+rxx;
end
Rx=Rxx/N;
%做特征值分解
[V,D]=eig(Rx);
p=4;
%对特征值和特征向量进行处理
for i=1:m;
    d(i)=D(m+1-i,m+1-i);
end
V1=zeros(m);
for i=1:m;
    V1(1:m,m+1-i)=V(1:m,i);
end
%取信号子空间
% l=m;
% for i=1:m-1;
%    dd=d(i)-d(i+1);
%    if dd>l,
%       l=i;
%        break
%    end
% end
s=zeros(m,p);
for i=1:p,
   s(1:m,i)=V1(1:m,i);
end
%用信号子空间公式求功率谱
df=0.001;
I=eye(m);
f=0.001:0.001:0.5;
for i=1:length(f);
    wi=2*pi*df*i;
    for t=1:m;
        aw(t)=exp(-j*(t-1)*wi);
    end
    p(i)=1/(aw*(I-s*s')*aw');
end
pw=10*log10(abs(p));
plot(f,pw);