gusucode.com > MVDR算法仿真源码程序 > MVDR算法仿真源码程序/code/MVDR_new.m
function [MSC]=MVDR_new(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);