doa_classical.txt

这是利用MUsic和盲分离算法求解DOA的程序

%%%%%%%%-2008.5.19-%%%%%%%%%
%%%%%%%%---classical DOA based on subband ICA---%%%%%%%%%
clc;
clear;
N=20000;
S=wavread('D:\MATLAB\R2006a\work\old work\新实验2008\wdrums_liverec_5cm_mix.wav');
X1=S(:,1)';X2=S(:,2)';
X1=X1(1:N);
X1=X1-mean(X1);
X2=X2(1:N);
X2=X2-mean(X2);     
X=[X1;X2]';
%-------------windowed Fourier transform-------------%
N2=512;      %频点数
OVERLAP=256;
N3=floor((N-N2)/OVERLAP)+1;       %窗移动数 N3路信号，每路N2点
fs=16000;
win=zeros(1,N2);
for k=1:N2
   win(k)=0.54-0.46*cos(2*pi*k/N2);%--汉明窗--%
   F(k)=fs*(k-1)/N2;
end                                      
meshw=zeros(1,N2);
for K=1:N2
   meshw(K)=K;
end
ts=zeros(1,N3);        %------ts记录每次窗移位置 即窗移动后 窗的起点是信号点的哪个点------%
xw1=zeros(N3,N2);
xw2=zeros(N3,N2);
for i=1:N3
    if i==1
        j=1;
    else
        j=(i-1)*OVERLAP+1;
    end
    ts(i)=j;
    z1=X1(j:j+N2-1);
    z2=X2(j:j+N2-1);
    xw1(i,:)=fft(z1.*win,N2);
    xw2(i,:)=fft(z2.*win,N2);
end

%----------separation on each frequency bin--------------
for n1=1:N2
    z(:,1)=xw1(:,n1);
    z(:,2)=xw2(:,n1);         
%-----whitten before MUSIC-----%
%     B1=conj(xw1(:,n1)');
%     B2=conj(xw2(:,n1)');
%     zz=[B1;B2];
%     Cx=cov(zz');
%     [e,d]=eig(Cx);
%     w0=inv(sqrt(d))*e';
%     Q{n1}=w0;
%     zz=w0*zz;
%     z(:,1)=zz(1,:);
%     z(:,2)=zz(2,:);
% %%%%%%-----music-----%%%%%%%%
  z  = z - ones(N3,1) * mean(z);
  R = conj(z' * z) / N3;
  [E,D]=eig(R);
  dspace=0.05;
  theta = [-90:1:90]' * (pi/180); 
  omega = 2*pi*sin(theta);
  mth   = length(theta); 
      for k=1:1:mth
          omega(k)=omega(k)*dspace*F(n1)/343;
          steer= exp(-sqrt(-1) * omega(k) * [0:1]');
          BP1(n1,k)=(abs(E(:,1)'*steer))^2;
          BP2(n1,k)=(abs(E(:,2)'*steer))^2;
      end 
%%%%%%%%%%%%%%%%%%    
    B1=conj(xw1(:,n1)');
    B2=conj(xw2(:,n1)');
    Z=[B1;B2];
    Cx=cov(Z');
    [e,d]=eig(Cx');
    w0=inv(sqrt(d))*e';
    Q{n1}=w0;
    Z=w0*Z;
    [W,y]=jade(Z,2);
    Y{n1}=y;                                           
    B{n1}=inv(W);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%---------------------%
for i=1:N2
    wf=B{i};%*inv(Q{i});
    w1=wf(1,:);%w1=w1/norm(w1);
	w2=wf(2,:);%w2=w2/norm(w2);
    
    dspace=0.05;
    theta = [-90:1:90]' * (pi/180); 
    omega = 2*pi*sin(theta);
    mth   = length(theta); 
    for k=1:1:mth
          omega(k)=omega(k)*dspace*F(i)/343;
          steer= exp(-sqrt(-1) * omega(k) * [0:1]');
          BP3(i,k)=abs(w1*steer);
          BP4(i,k)=abs(w2*steer);
    end 
end


figure(1);
plot(theta,20*log10(BP3(1,:)),'r');
hold on
plot(theta,20*log10(BP4(1,:)),'g');

for i=20:125
   plot(theta,20*log10(BP3(i,:)),'r');
   hold on
   plot(theta,20*log10(BP4(i,:)),'g');
   hold on
end