QQ登录

只需要一步,快速开始

 注册地址  找回密码
查看: 4501|回复: 0
打印 上一主题 下一主题

麦克风阵列声源定位 SRP-PHAT(二)

[复制链接]
字体大小: 正常 放大
浅夏110 实名认证       

542

主题

15

听众

1万

积分

  • TA的每日心情
    开心
    2020-11-14 17:15
  • 签到天数: 74 天

    [LV.6]常住居民II

    邮箱绑定达人

    群组2019美赛冲刺课程

    群组站长地区赛培训

    群组2019考研数学 桃子老师

    群组2018教师培训(呼伦贝

    群组2019考研数学 站长系列

    跳转到指定楼层
    1#
    发表于 2020-5-15 15:04 |只看该作者 |倒序浏览
    |招呼Ta 关注Ta |邮箱已经成功绑定
    DOA
    * ^. Y! o9 Y7 ^  声源定位方法一般可分为三类,一种是基于TDOA的两步算法(two-stage algorithm),一种是基于空间谱估计如MUSIC等,还有就是基于beamforming的方法,也就是这里要介绍的可控波束响应(steered-response power),
      k* s2 T7 k8 P- N# }8 h
    . B( c% E4 W# R$ f; v. N. O- W* jsteered-response power. k  S6 h8 @. d
      可控波束响应是利用波束形成(beamforming)的方法,对空间不同方向的声音进行增强,得到声音信号最强的方向就被认为是声源的方向。
    , t% b- C, c: B* j. J  B' i  上一篇中简单介绍了麦克风阵列的背景知识,最简单的SRP就是利用延时-累加(delay-and-sum)的方法,寻找输出能量最大的方向。 3 i% Q' p" Q+ R- e. C4 ]
      其中,语音信号为宽带信号,因此需要做宽带波束形成,这里我们在频域实现
    ) H% d1 }( [3 K1 b7 O
    0 x$ Q( K, k# M- e频域宽带波束形成
    + d, h+ i4 p+ Z, k. \) j/ o  频域宽带波束形成可以归类为DFT波束形成器,结构如下图
    ' P! d5 [( j3 i$ t  `4 ~0 Y& c0 ^& F2 H9 f- D
    . R. P+ V7 e& R$ a# O$ ^

    频域处理也可以看做是子带处理(subband),DFT和IDFT的系数分别对应子带处理中的分析综合滤波器组,关于这一种解释,可参考《传感器阵列波束优化设计与应用》第六章。

    频域宽带延时累加波束形成的基本过程就是信号分帧加窗->DFT->各频点相位补偿->IDFT
    : Q4 Y- [- I% l0 u代码实现如下


    2 ~- U+ k3 x0 `nction [ DS, x1] = DelaySumURA( x,fs,N,frameLength,inc,r,angle)
    : n# P0 ]' m6 L. o) q: m/ K%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5 g4 R3 }- j; a1 x
    %frequency-domain delay-sum beamformer using circular array1 J8 z! R8 F. i4 b5 I+ c; p0 j9 N
    %   5 E9 B1 H; Z: f( m, s" W
    %      input :
    + l/ W; O' D, [+ {+ A, r% {0 [& n%          x : input signal ,samples * channel. }1 R" r5 I  t+ v
    %          fs: sample rate
    ! K, b8 ^9 p3 V. R%          N : fft length,frequency bin number/ l0 B% f: r! k8 f5 D
    %frameLength : frame length,usually same as N
    7 k6 z) K# t: S& S0 {" r4 y% N%        inc : step increment
    ! _. w9 ]' b4 [4 q& _% a%          r : array element radius
    . J* j$ L; l  i. r$ C$ K" G) j%      angle : incident angle: j) V3 k/ a/ O8 o8 ]# Y
    %
    7 O0 O8 ]5 S9 u  ~) g%     output :
    . A! F# [6 |8 I6 H) d6 f( ^' b%         DS : delay-sum output- K8 V) \" o0 F" p1 @3 A, B3 Q
    %         x1 : presteered signal,same size as x
    ; \1 {6 ?. Z  p9 K8 O* |9 N%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    7 C$ T) S6 I3 O( z3 P3 ]* T3 m6 J/ m/ [# R% X0 M
    c = 340;, }* W0 x' z6 w' F
    Nele = size(x,2);7 [2 R# ?0 a3 l+ B: |4 S" G& ^8 l
    omega = zeros(frameLength,1);8 }/ ^6 s* K' k% l! i
    H = ones(N/2+1,Nele);
    # m% @% |+ N/ F2 [
    4 P4 g$ U1 q% P& ^& ~) }! d+ ytheta = 90*pi/180; %固定一个俯仰角3 W( K8 g( ]! ]* J* v7 V' \  G
    gamma = [30 90 150 210 270 330]*pi/180;%麦克风位置
    $ H9 U& z8 Y4 h; Xtao = r*sin(theta)*cos(angle(1)-gamma)/c;     %方位角 0 < angle <360# i" B5 W. |3 l# O& x' P
    yds = zeros(length(x(:,1)),1);% ]/ x7 m7 j+ O8 L' J- x
    x1 = zeros(size(x));, |3 U  b/ x! G, S

    ! M  W, O* m, k, c/ g* p% frequency bin weights% E" V7 E* [3 F0 B% q
    % for k = 2:1:N/2+15 }8 k1 k( i' T( x* u4 t
    for k = 1:1:5000*N/fs7 h* z; p6 S9 u1 W+ |& C
        omega(k) = 2*pi*(k-1)*fs/N;   
    6 P1 v- L& u% p( T4 E. T- g$ f$ D    % steering vector! o/ p% x$ [; K9 X: A" T
        H(k, = exp(-1j*omega(k)*tao);
    ) W/ H/ T3 d+ x& o3 r# n+ f& S3 {3 k1 `' mend' T/ u& {+ }3 L% o1 a: o
    , V* F+ _# n/ T
    for i = 1:inc:length(x(:,1))-frameLength
    % Z- m9 R! \$ Y, d( p$ O9 T# x' B, d) A1 ]( s3 y) a
        d = fft(bsxfun(@times, x(i:i+frameLength-1,,hamming(frameLength)));
    ! l2 K% e6 b# Q' m" d' x' s* w- D: m% |% h
        x_fft=bsxfun(@times, d(1:N/2+1,,H);+ g/ k# O. q7 e! X0 \* C( v
    8 D4 c* i  a+ P
        % phase transformed8 `" }$ }& J# c- B3 Z% p
        %x_fft = bsxfun(@rdivide, x_fft,abs(d(1:N/2+1,));
    : ?9 h$ t! e- `$ B0 F. x    yf = sum(x_fft,2);
    + t4 Y$ U7 {+ d7 L9 m5 w; x    Cf = [yf;conj(flipud(yf(2:N/2)))];' D  }& t; l1 h+ ?$ m4 p

    * `& T0 I. k# }4 f; o1 F4 `    % 恢复延时累加的信号2 K4 I! S. R( [  Q0 @
        yds(i:i+frameLength-1) = yds(i:i+frameLength-1)+(ifft(Cf));4 q# n% ]# @, D- x( k
    - P: z1 C2 b6 {' H$ w
        % 恢复各路对齐后的信号7 ]. i8 U# ]" k' }+ Y4 L" m4 @
        xf  = [x_fft;conj(flipud(x_fft(2:N/2,))];
    & P1 ]' g" E( c$ e    x1(i:i+frameLength-1, = x1(i:i+frameLength-1,+(ifft(xf));
    $ H/ B* \! x) {5 h4 E9 A. }! i  f- Wend
    / ?5 p- C( |5 s$ n; A7 eDS = yds/Nele;  
    0 w+ N( k, E: {  y8 h
    " I) K; T# z6 Z* Z! x( [' t8 fend
    4 i$ L. K4 J9 O2 k- N& @9 g然后遍历各个角度重复调用这个函数,测试实际录音数据,代码如下( b) [6 |! D9 S4 m" {
    6 j- o% z! J* A
    %% SRP Estimate of Direction of Arrival at Microphone Array
      u& }3 C0 X  G% Frequency-domain delay-and-sum test6 P  |% G  j  o3 |
    %  - G* ~- h2 ~5 Z0 k) L
    %%
      [4 D/ a( G3 S* {/ s1 u$ Z% @
      M' B) K. r7 A: q0 O& ?$ F4 u% x = filter(Num,1,x0);( s( P/ K  C# T
    c = 340.0;4 h0 i* [3 c- E* q; l8 K* s5 R
    ! x: w: ^9 A; p5 D
    % XMOS circular microphone array radius) ?& x  ^: S# C7 z0 g
    d = 0.0420;
    / R" o) c2 M( I8 a6 n5 Y$ S% more test audio file in ../../TestAudio/ folder
    - s' l! Z8 _$ }+ ypath = '../../TestAudio/XMOS/room_mic5-2/';
    , k  e. I' z2 H: Q  E$ G5 C* ~[s1,fs] = audioread([path,'音轨-2.wav']);
    , a- z/ O1 w' s6 ^0 y) w8 Os2 = audioread([path,'音轨-3.wav']);6 O4 ]0 p( r# F! i* D
    s3 = audioread([path,'音轨-4.wav']);$ B: Q" J! V: B1 A' l
    s4 = audioread([path,'音轨-5.wav']);7 K* t. N$ t: B* I& `
    s5 = audioread([path,'音轨-6.wav']);
    1 ]2 i7 E* l; }4 f% d) O1 G9 X, Js6 = audioread([path,'音轨-7.wav']);
    1 k& T. Y& j- x: V( I  rsignal = [s1,s2,s3,s4,s5,s6];
      z  c" B: G) O; |5 \4 ]M = size(signal,2);5 S& _3 f! u( N
    %%
    + D& r, X+ i; ]5 i1 u, Jt = 0;  Y8 T: C1 F3 @. b9 Y5 f
    # Q7 ~. d2 c" _; E9 U
    % minimal searching grid
    - ~9 b' @$ \, E6 F5 m/ H. _/ @step = 1;
    5 x2 F% x, }* p- Y& Z4 I7 ]0 t; M. z* c2 [' s
    P = zeros(1,length(0:step:360-step));* n8 b6 v8 l: Q+ f& _' ]# v
    tic
    ' N$ r0 Y  X" x3 \4 Bh = waitbar(0,'Please wait...');# n% x6 |8 t. J$ a
    for i = 0:step:360-step* R/ X, w* `0 d* @3 x6 q1 `
        % Delay-and-sum beamforming. ], {0 q# z0 H! \* ~
        [ DS, x1] = DelaySumURA(signal,fs,512,512,256,d,i/180*pi);
    2 L# f+ F8 O  T4 P    t = t+1;
    ) i: }* Q5 N/ ^; Z    %beamformed output energy
    3 \  z! P& j! m% |6 d8 d; h    P(t) = DS'*DS;' Y, Y* r! f$ O3 b. I
        waitbar(i / length(step:360-step))
    7 r- K% Q3 @. P) f5 fend4 U9 u/ l4 n7 Q* y8 L6 A
    toc
    ; D3 W0 L1 y  {& bclose(h)   J9 x" ^3 \/ H, D
    [m,index] = max(P);6 U6 ^# r* l/ E% J3 x- }' ^
    figure,plot(0:step:360-step,P/max(P))4 Y* u7 \% ^- {& \$ _, S0 f
    ang = (index)*step
    ) x+ V* V5 a$ o# e" m& w1 w2 r$ T' J+ `" i5 I& w
    程序中用的是圆阵,可以进行二维方向角扫描,不过这里为了简便就固定了俯仰角,只扫描方位角,结果如下
    : {3 o7 H5 J/ q9 k6 g
    5 l1 A) H; }4 ]9 Z$ K* Q, T

    结果与预期相同

    PHAT加权

      与GCC-PHAT方法相同,这里也可以对幅度做归一化,只保留相位信息,使得到的峰值更明显,提高在噪声及混响环境下的性能 # Z& [9 n8 j, K/ w
      上面代码中加上这一句

    * M5 q! f$ b7 g' \. u
    %x_fft = bsxfun(@rdivide, x_fft,abs(d(1:N/2+1,));
    9 j, B/ H% I# R4 m测试同样的文件,结果如下 9 b; Q2 A; ?# z5 A. G

    $ d4 c$ m6 m+ ]% I* A5 ~; b. l对比可以看到,PHAT加权的方法性能更好" r- e4 T: H3 \
    参考! W1 Z' G7 Z% y
    1.《SRP-PHAT-A High-Accuracy, Low-Latency Technique for Talker Localization in Reverberant Environments Using Microphone Arrays》 ) d  S* w% l! n8 N/ `* \5 V6 e
    2. 《传感器阵列波束优化设计与应用》
    + }2 q: M# |% t————————————————$ k2 m# K: q" L- _/ }
    版权声明:本文为CSDN博主「373955482」的原创文章,遵循CC 4.0 BY-SA版权协议,转载请附上原文出处链接及本声明。: z3 c* J/ @9 E, b3 j
    原文链接:https://blog.csdn.net/u010592995/article/details/81586504
    ! u' n* `, O9 C! ^' x
    : m$ w5 E  h* x- l: M) C+ e1 e! C: C

    1 i8 A# r7 M) C& b) x
    , m1 E' j* K& V  F& b
    ! M% U0 L! D: z' k0 i! o4 @
    zan
    转播转播0 分享淘帖0 分享分享0 收藏收藏0 支持支持0 反对反对0 微信微信
    您需要登录后才可以回帖 登录 | 注册地址

    qq
    收缩
    • 电话咨询

    • 04714969085
    fastpost

    关于我们| 联系我们| 诚征英才| 对外合作| 产品服务| QQ

    手机版|Archiver| |繁體中文 手机客户端  

    蒙公网安备 15010502000194号

    Powered by Discuz! X2.5   © 2001-2013 数学建模网-数学中国 ( 蒙ICP备14002410号-3 蒙BBS备-0002号 )     论坛法律顾问:王兆丰

    GMT+8, 2026-4-21 07:40 , Processed in 0.304088 second(s), 51 queries .

    回顶部