add test code

Author: wangwei2009

README.md

@@ -1,2 +1,2 @@
-# GCC-PHAT
+# DOA
 Estimate the direction of arrival using the SRP-PHAT and GCC-PHAT algorithm.

matlab/ArrayTest/DS_ULA.m

@@ -0,0 +1,18 @@
+M=8;                 % 阵元数
+c=340;                         % 声速
+f = 2125;                      % 信号频率
+d = 0.08;                      % 阵元间距
+theta=linspace(-pi/2,pi/2,200);% 入射信号角度范围
+theta0=30*pi/180;                      % 注视方向
+w=exp(1j*2*pi*f*sin(theta0)*[0:M-1]'*d/c);       % 导向向量
+p = zeros(length(1:length(theta)),1);
+for  j=1:length(theta)                                     % 角度扫描
+    a=exp(-1j*2*pi*f*sin(theta(j))*[0:M-1]'*d/c);% 入射信号方向向量
+    p(j)=sum(w.*a)/M;                                        % 延时-求和
+end
+% 画图
+figure;
+plot(theta/pi*180,abs(p)),grid on
+xlabel('degree')
+ylabel('amplitude')
+title('8阵元均匀线阵方向图')

matlab/ArrayTest/DS_URA_2D.m

@@ -0,0 +1,26 @@
+M=8;                 % 阵元数
+c=340;                         % 声速
+f = 2125;                      % 信号频率
+r = 0.08;                      % 阵元间距
+phi=linspace(0,2*pi,360);% 入射信号角度范围
+
+theta = 80*pi/180; %固定一个俯仰角
+gamma = (0:360/M:360-360/M)*pi/180;%麦克风位置
+
+tao = r*sin(theta)*cos(angle(1)-gamma)/c;     %方位角 0 < angle <360
+
+phi0=60*pi/180;                      % 注视方向
+w=exp(1j*2*pi*f*r*cos(phi0-gamma)*sin(theta)/c);       % 导向向量
+p = zeros(length(1:length(phi)),1);
+for  j=1:length(phi)                                     % 角度扫描
+    a=exp(-1j*2*pi*f*r*cos(phi(j)-gamma)*sin(theta)/c);% 入射信号方向向量
+    p(j)=sum(w.*a)/M;                                        % 延时-求和
+end
+% 画图
+figure;
+plot(phi/pi*180,abs(p)),grid on
+xlabel('degree')
+ylabel('amplitude')
+title('均匀圆阵方向图')
+
+figure,polarplot(phi,abs(p))

matlab/ArrayTest/DS_URA_3D.m

@@ -0,0 +1,29 @@
+M=8;                 % 阵元数
+c=340;                         % 声速
+f = 2125;                      % 信号频率
+r = 0.08;                      % 阵元间距
+phi=linspace(0,2*pi,360);% 入射信号方位角范围
+
+theta = linspace(0,pi/2,90);% 入射信号俯仰角范围
+
+gamma = (0:360/M:360-360/M)*pi/180;%麦克风位置
+
+
+phi0=130*pi/180;                      % 注视方向
+theta0 = 80*pi/180;
+
+w=exp(1j*2*pi*f*r*cos(phi0-gamma)*sin(theta0)/c);       % 导向向量
+p = zeros(length(1:length(phi)),length(theta));
+for i = 1:length(theta)
+    for  j=1:length(phi)                                     % 角度扫描
+        a=exp(-1j*2*pi*f*r*cos(phi(j)-gamma)*sin(theta(i))/c);% 入射信号方向向量
+        p(j,i)=sum(a.*w)/M;                                        % 延时-求和
+    end
+end
+% 画图
+% figure;
+% plot(phi/pi*180,abs(p)),grid on
+% xlabel('degree')
+% ylabel('amplitude')
+% title('均匀圆阵方向图')
+figure,mesh(abs(p))

matlab/README.md

@@ -0,0 +1,8 @@
+# A simple DOA GUI
+a simple GUI to show the direction of arrival of a signal
+
+
+## usage
+1. run DispAngle.m
+2. select your multi-channel soundcard number in command window
+3. click Update button to see the real time direction

matlab/SRP-PHAT/DelaySumURA.m

@@ -19,21 +19,21 @@
 c = 340;
 Nele = size(x,2);
 omega = zeros(frameLength,1);
-H = zeros(N/2+1,Nele);
+H = ones(N/2+1,Nele);
 
 theta = 85*pi/180; %固定一个俯仰角
 gamma = [0 90 180 270]*pi/180;%麦克风位置
 tao = r*sin(theta)*cos(angle(1)-gamma)/c;     %方位角 0 < angle <360
 yds = zeros(length(x(:,1)),1);
 x1 = zeros(size(x));
 for i = 1:inc:length(x(:,1))-frameLength
-    for k = 2:N/2+1
+    for k = 2:4:N/2+1
         omega(k) = 2*pi*(k-1)*fs/N;   
 
         % steering vector
         H(k,:) = exp(-1j*omega(k)*tao);
     end
-    d = fft(bsxfun(@times, x(i:i+frameLength-1,:),hamming(frameLength)));
+    d = fft(bsxfun(@times, x(i:i+frameLength-1,:),hamming(frameLength)'));
 %     d = fft(x(i:i+frameLength-1,:).*hamming(frameLength)');
 %     x_fft = d(1:129,:).*H;%./abs(d(1:129,:));
     x_fft=bsxfun(@times, d(1:N/2+1,:),H);