mvdr.m

来自「窄带信号mvdr波束形成,这几天的学习结果」· M 代码 · 共 33 行

clear all;
cj=sqrt(-1);

N=256;
dt=0.001;
df=1/(N*dt);
f1=200;
f2=350;
f3=100;
t=dt*(-N/2:N/2-1);
f=df*(-N/2:N/2-1);

% Use Fourier method
s=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t)+0.8*sin(2*pi*f3*t+pi/3);
s=s+randn(1,N)*0.01;
S=fftshift(fft(fftshift(s)));
figure; plot(f,abs(S));

% Use MVDR method
R=zeros(N,N);
for i=1:3000
s=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t)+0.8*sin(2*pi*f3*t+pi/3);
s=s+randn(1,N)*0.01;
R=R+s'*s;
end
R=R/3000;
Rinv=inv(R);
P=zeros(1,N);
for i=1:N
s=exp(cj*(0:N-1)*2*pi*(f(i)*dt)).';
P(i)=1/abs(s'*Rinv*s);
end
figure; plot(f,P);