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:3for k=1:gun=awgn(Un,snr(i),'measured');%% SVDA=zeros(M,N-M+1);       %构造数据矩阵for n=1:N-M+1A(:,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:Pfor m=1:Ma(m,kk)=exp(-1j*omega(kk)*(m-1));endendPmvdr=zeros(1,P);       %MVDR谱for kk=1:PPmvdr(kk)=1/(a(:,kk)'*invphi*a(:,kk));endPmvdr=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