clear ;close all;clc;

产生输入信号

N = 1024;               %样本点数
snr=[20 25 30];         %信噪比
n=0:N-1;                %数据轴
g=100;                  %蒙特卡诺仿真次数
M=14;                   %阶数
Pmvdr_s=zeros(3,1024);  %存放MVDR谱
signal1 = exp(1i*0.15*2*pi*n + 1i*2*pi*rand);
signal2 = exp(1i*0.25*2*pi*n + 1i*2*pi*rand);
signal3 = exp(1i*0.30*2*pi*n + 1i*2*pi*rand);
Un = signal1 + signal2 + signal3;

蒙特卡诺仿真

for i=1:3
    for k=1:g
        un=awgn(Un,snr(i),'measured');
        %% SVD
        A=zeros(M,N-M+1);       %构造数据矩阵
        for n=1:N-M+1
            A(:,n)=un(M+n-1:-1:n);
        end
        [U,S,V]=svd(A');
        invphi=V*inv(S'*S)*V';  %类似于自相关矩阵逆矩阵
        %% 计算MVDR谱
        P = 1024;               %MVDR扫描点数
        f=linspace(-0.5,0.5,P); %频率轴
        omega=2*pi*f;           %归一化角频率
        a=zeros(M,P);           %频率导向矢量
        for kk=1:P
            for m=1:M
                a(m,kk)=exp(-1j*omega(kk)*(m-1));
            end
        end
        Pmvdr=zeros(1,P);       %MVDR谱
        for kk=1:P
            Pmvdr(kk)=1/(a(:,kk)'*invphi*a(:,kk));
        end
        Pmvdr=10*log10(abs(Pmvdr/max(abs(Pmvdr))));
        Pmvdr_s(i,:) =Pmvdr_s(i,:)+Pmvdr;
    end
end
Pmvdr_s=Pmvdr_s/g;

绘图

figure;
hold on
plot(f,Pmvdr_s(1,:));
plot(f,Pmvdr_s(2,:));
plot(f,Pmvdr_s(3,:));
hold off
xlabel('w/2Π');ylabel('归一化功率谱 (dB)');
legend('SNR=20','SNR=25','SNR=30');
title('MVDRSVD频率估计方法');
grid on