【雷达】阵列信号处理MVDR程序高分辨方位估计附matlab代码-CFANZ编程社区

1 简介

The minimum variance distortionless response (MVDR),originally developed by Capon for frequency-wavenumberanalysis, is a very well established method in array processing. It is also used in spectral estimation. The aim of thispaper is to show how the MVDR method can be used to estimate the magnitude squared coherence (MSC) function,which is very useful in so many applications but so fewmethods exist to estimate it. Simulations show that our algorithm gives much more reliable results than the one basedon the popular Welch’s method.

2 完整代码

function [MSC]=coherence_MVDR(x1,x2,L,K);
%% This program computes the coherence function between 2 signals 
%% x1 and x2 with the MVDR method.
%% This algorithm is based on the paper by the same authors:
%% J. Benesty, J. Chen, and Y. Huang, "A generalized MVDR spectrum," 
%% IEEE Signal Processing letters, vol. 12, pp. 827-830, Dec. 2005.
%% x1, first signal vector of length n
%% x2, second signal vector of length n
%% L is the length of MVDR filter or window length
%% K is the resolution (the higher K, the better the resolution)
%initialization
xx1     = zeros(L,1);
xx2     = zeros(L,1);
r11     = zeros(L,1);
r22     = zeros(L,1);
r12     = zeros(L,1);
r21     = zeros(L,1);
%construction of the Fourier Matrix
F       = zeros(L,K);
l       = [0:L-1]';
f       = exp(2*pi*l*j/K);
for k = 0:K-1
    F(:,k+1) = f.^k;
end
F       = F/sqrt(L);
%number of samples, equal to the lenght of x1 and x2
n       = length(x1);
for i = 1:n
    xx1 = [x1(i);xx1(1:L-1)];
    xx2 = [x2(i);xx2(1:L-1)];
    r11 = r11 + xx1*conj(xx1(1));
    r22 = r22 + xx2*conj(xx2(1));
    r12 = r12 + xx1*conj(xx2(1));
    r21 = r21 + xx2*conj(xx1(1));
end
r11 = r11/n;
r22 = r22/n;
r12 = r12/n;
r21 = r21/n;
%
R11 = toeplitz(r11);
R22 = toeplitz(r22);
R12 = toeplitz(r12,conj(r21));
%
%for regularization
Dt1     = 0.01*r11(1)*diag(diag(ones(L)));
Dt2     = 0.01*r22(1)*diag(diag(ones(L)));
%
Ri11    = inv(R11 + Dt1);
Ri22    = inv(R22 + Dt2);
Rn12    = Ri11*R12*Ri22;
%
Si11    = real(diag(F'*Ri11*F));
Si22    = real(diag(F'*Ri22*F));
S12     = diag(F'*Rn12*F);
%
%Magnitude squared coherence function
MSC     = real(S12.*conj(S12))./(Si11.*Si22);

clear all
%% This program calls the main program "coherence_MVDR",
%% which computes the coherence function with the MVDR method
%example on how the coherence function works with the MVDR method
%see papers for more explanations
n       = 1024;         %number of samples
nT      = [0:n-1]';     %time axis
Nf      = 5;            %number of frequencies of high coherence
f       = zeros(Nf,1);
f(1) = 0.05; f(2) = 0.06; f(3) = 0.07; f(4) = 0.08; f(5) = 0.09;
%
fw      = 2*pi*f;
x1      = randn(n,1); %first signal
x2      = randn(n,1); %second signal
for i = 1:Nf
    x1  = x1 + cos(fw(i)*nT);
    x2  = x2 + cos(fw(i)*nT + 2*pi*rand(1,1));
end
%coherence with MATLAB
WL      = 100; %window length
[cx1x2,w] = cohere(x1,x2,WL);
%
figure(1);
%
subplot(2,1,1)
plot(w/2,cx1x2)
%legend('');
grid on;
ylabel('MSC');
title('(a)')
axis([0 1/2 0 1]);
%Coherence with the MVDR method
K = 200; %to increase resolution
L = 100; %window length
[MSC]=coherence_MVDR(x1,x2,L,K);
%
K2 = K/2;
MSCf = MSC(1:K2);
wwf = [0:1/K2:1-1/K2]';
%
subplot(2,1,2)
plot(wwf/2,MSCf)
grid on;
ylabel('MSC');
xlabel('Frequency');
title('(b)')
axis([0 1/2 0 1]);