music.m

来自「ar模型 bt算法 ls rls lms music等数字信号处理的源代码」· M 代码 · 共 60 行

clc
clear
close
n=1:1000;
w0=1.3*pi;
x=exp(j*pi*n-j*pi)+exp(j*w0*n-j*.6*pi)+sqrt(.1)*randn(1,length(n));

%%%%%%%%%%%%%   Estimate auto-realtion   %%%%%%%%%%%%%

r=zeros(1,501);
for m=1:501
    for n=m:1000
        r(m)=r(m)+x(n)*conj(x(n-m+1));
    end
end
r=r/1000;

%%%%%%%%%%%   Estimate classical  PSD    %%%%%%%%%%%%%

r1=[conj(fliplr(r(2:501))) r];
k=0:800;
n1=-500:500;
P_BT=r1*exp(-j*pi/400).^(n1'*k);
figure
plot(pi/400*k,10*log10(P_BT/max(P_BT)))
xlabel('w')
ylabel('normalized PSD(db)')

%%%%%%%% construct 4th order correlation matrix      %%%%%%%%

R=zeros(4);
for m=1:4
    for n=1:4
        if(m<=n) R(m,n)=r(n-m+1);
        else R(m,n)=conj(r(m-n+1));
        end
    end
end
%%%%%%%%%%   music algorithm     %%%%%%%%%%%%%%%%%%%%%%%%

[V D]=eig(R);
w=linspace(0,2*pi,300)';
a=exp(-j*[0 1 2 3]'*w');
pmusic=1./((abs(a'*V(:,1))).^2+(abs(a'*V(:,2))).^2);
figure
plot(w,10*log10(pmusic))
xlabel('w')
ylabel('P_MUSIC(db)')

%%%%%%%%%%%%%%%%%%%%%%  MVDR    %%%%%%%%%%%%%%%%

for m=1:300
    a1=[1 exp(-j*w(m)) exp(-2*j*w(m)) exp(-3*j*w(m)) ].';
    pmvdr(m)=1./(a1'*inv(R)*a1);
end
figure
plot(w,10*log10(pmvdr))
xlabel('w')
ylabel('P_MVDR(db)')