QQ登录

只需要一步,快速开始

 注册地址  找回密码
查看: 4505|回复: 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 |邮箱已经成功绑定
    DOA1 L3 }$ c+ R3 l3 {& X- o- H) {
      声源定位方法一般可分为三类,一种是基于TDOA的两步算法(two-stage algorithm),一种是基于空间谱估计如MUSIC等,还有就是基于beamforming的方法,也就是这里要介绍的可控波束响应(steered-response power),- o3 \5 e# y. i, u/ S. ?0 h

    ' O* k7 m9 c, u: s; s8 f' bsteered-response power3 H; q  `" ?$ r+ S1 n
      可控波束响应是利用波束形成(beamforming)的方法,对空间不同方向的声音进行增强,得到声音信号最强的方向就被认为是声源的方向。
      s4 l6 f, z+ Y) W3 g1 w! q' e$ V  上一篇中简单介绍了麦克风阵列的背景知识,最简单的SRP就是利用延时-累加(delay-and-sum)的方法,寻找输出能量最大的方向。
    2 T' L% u0 {1 i3 l  其中,语音信号为宽带信号,因此需要做宽带波束形成,这里我们在频域实现
    6 {+ \. ^" C& v
    . I  A# T& d8 {6 `频域宽带波束形成
    : `! U' R8 L" F# P  频域宽带波束形成可以归类为DFT波束形成器,结构如下图 / d- t& W( d2 z8 m  J/ ~) R; T5 E

    # n  v6 D: F  U* n# O. ?. [0 U9 _# W7 g0 ^# F

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

    频域宽带延时累加波束形成的基本过程就是信号分帧加窗->DFT->各频点相位补偿->IDFT 9 F+ t7 M7 l' v, [0 T8 n7 b
    代码实现如下


    * S* D0 B0 v* R3 o: O2 enction [ DS, x1] = DelaySumURA( x,fs,N,frameLength,inc,r,angle)
    $ b" b5 e. _; _# D( T7 F0 G0 u%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % {& w* Z7 ]) T. O9 ?%frequency-domain delay-sum beamformer using circular array; Q3 [, K/ {# V5 m& Y* l3 O
    %   
    7 x& k% J& K2 ~%      input :
    # S" L3 ~4 D4 d. z5 A9 d& s%          x : input signal ,samples * channel8 i, D* |; S& D8 h! q. a: P. H
    %          fs: sample rate
    $ F- K# o) d. E- F9 v; M& ]' A: L%          N : fft length,frequency bin number) S, d$ r, g+ Z6 H# ^, {9 e2 c
    %frameLength : frame length,usually same as N
    ( T% z; U6 c* Z( s$ i; e%        inc : step increment/ q' f( B- d2 a. _% k
    %          r : array element radius. d4 Q9 a; Y- ]/ t! f
    %      angle : incident angle
    5 Z  P1 ~- e0 t$ H; U%
    9 l  v2 Z" f3 I1 ^8 ]# i) i%     output :
    " c5 }. A; T, \. C( ^%         DS : delay-sum output
    ) T- m% i3 k  E%         x1 : presteered signal,same size as x  y/ K  G! j, R4 V# I( H/ W% X
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    4 W: q/ E! _  k, y- y
    ; j- n4 H3 {$ j. q* nc = 340;
    $ |7 {- C! |* n1 TNele = size(x,2);2 l7 ^( P' q" F! g* w. J
    omega = zeros(frameLength,1);. f% C7 Z) Z1 G1 G8 W3 a$ }# `" h
    H = ones(N/2+1,Nele);
    + x6 Z  {  G6 M8 T
    6 F9 a5 x* f1 W" U) b7 T# itheta = 90*pi/180; %固定一个俯仰角8 ]  ^6 e6 [. X/ c# L% l' G
    gamma = [30 90 150 210 270 330]*pi/180;%麦克风位置
    5 j/ o- c5 K4 b, v* o; A9 M  xtao = r*sin(theta)*cos(angle(1)-gamma)/c;     %方位角 0 < angle <360
    6 q, {. T7 W( y( \4 X9 Oyds = zeros(length(x(:,1)),1);
    - t$ F5 N- i2 X1 E, y/ Cx1 = zeros(size(x));
    : N( w, K5 V% N- O) Y' w* r% a
    ) Y9 H" V: n$ e' f3 H" B% frequency bin weights2 q) {8 L, f* a" j; b' Z. H
    % for k = 2:1:N/2+14 w' [) a- m0 N* l. h! e# s, i
    for k = 1:1:5000*N/fs2 F7 `% S2 d! l$ ^
        omega(k) = 2*pi*(k-1)*fs/N;   
    3 i; c' N  k: b, M' Q    % steering vector: s5 U2 w! q2 Z+ m0 P* M- C
        H(k, = exp(-1j*omega(k)*tao);8 Q# _! D3 A! d2 J* Y
    end
    * A7 G: `; A9 j0 C1 J7 |$ ^  ]) g6 t! O7 P4 M! D& [. N7 j2 L
    for i = 1:inc:length(x(:,1))-frameLength6 f( ?% U, ?! C

    ; B' B! ^, Q7 h5 I$ B" U: q    d = fft(bsxfun(@times, x(i:i+frameLength-1,,hamming(frameLength)));! F- B* o1 J& F. h/ F

    1 g. ~. Y) B% p! d% A2 q/ W    x_fft=bsxfun(@times, d(1:N/2+1,,H);
    # v: `- |, h3 w  d" `3 `6 E( h" E" W1 R
        % phase transformed2 k% r; u) C! x
        %x_fft = bsxfun(@rdivide, x_fft,abs(d(1:N/2+1,));/ ?, y  e/ a7 F$ h2 `. X5 I9 ^4 x
        yf = sum(x_fft,2);
    . W7 a3 d9 M+ q    Cf = [yf;conj(flipud(yf(2:N/2)))];; r/ i' `( Y) u3 h  E& \6 @# m& F
    * q: N- b' W6 ?2 G$ Y
        % 恢复延时累加的信号1 q1 W: x% M5 U1 [' R
        yds(i:i+frameLength-1) = yds(i:i+frameLength-1)+(ifft(Cf));
    2 p8 f: o/ W  |! z' E
    ! j+ \0 ~* q6 o    % 恢复各路对齐后的信号+ Z" x* E5 S1 b
        xf  = [x_fft;conj(flipud(x_fft(2:N/2,))];' |5 c+ Y  X) B3 B
        x1(i:i+frameLength-1, = x1(i:i+frameLength-1,+(ifft(xf));. j0 v( Y9 w( ?. G
    end
    ) e0 o, Q2 T$ }! g9 P% IDS = yds/Nele;  " ^* e5 V* b5 q% f% b

    9 M  T# w4 D. P+ E2 D/ [" wend
    # S) ^' P: m* K. p/ v. Q/ X然后遍历各个角度重复调用这个函数,测试实际录音数据,代码如下
    # X/ E1 e& K! ~
      `; O0 V4 Y- y  ~4 U9 d/ w6 ]%% SRP Estimate of Direction of Arrival at Microphone Array$ h( d+ y3 s, \& A& l
    % Frequency-domain delay-and-sum test
    4 N  T9 L0 P7 a+ u2 n%  ! f4 @1 c1 T6 C0 U4 z
    %%+ k. o8 r  f( M

    1 V) f. s' g, I8 U. p. }& P% x = filter(Num,1,x0);! i, j( o2 `0 [) ~+ g, t5 L+ \( Z
    c = 340.0;. S# f5 n, g0 v; h3 [

    - L) ^2 X( E3 ]9 h; r% XMOS circular microphone array radius& L0 H% S1 O7 u! u0 y
    d = 0.0420;8 V  W( O% o4 k0 ]/ Z2 k
    % more test audio file in ../../TestAudio/ folder
    8 |( c8 Y5 r& e  `5 I5 M( rpath = '../../TestAudio/XMOS/room_mic5-2/';& y2 M3 Z* N3 U0 C: m
    [s1,fs] = audioread([path,'音轨-2.wav']);
    , V2 E) X) U- Ms2 = audioread([path,'音轨-3.wav']);* C( v; F* x$ y3 l
    s3 = audioread([path,'音轨-4.wav']);
    & A  X4 }* p- K1 C8 K" ]7 ts4 = audioread([path,'音轨-5.wav']);
    & k0 [3 z; Z) H3 Y" ?% x* Us5 = audioread([path,'音轨-6.wav']);4 G& k! M% C5 {2 Y- z- |: I
    s6 = audioread([path,'音轨-7.wav']);
    3 q5 {0 K% g. O: M: V% qsignal = [s1,s2,s3,s4,s5,s6];# j% F% A1 N1 f$ q3 q" v/ l
    M = size(signal,2);
    3 ~* o# W& `1 o0 ^6 B%%  r; X; D) M5 D9 Q: }6 n$ E
    t = 0;; p4 i  J6 ^# n, t; n

    # n2 [% U) W& v& F+ [1 G# H% v% O" W% minimal searching grid# C2 l- ~0 Y- ~+ X$ f6 }2 V
    step = 1;
    ( q. a4 y/ ^. {+ f7 C( ]* V  D2 R) j6 [  S
    P = zeros(1,length(0:step:360-step));! [9 L: _( i2 u$ Y2 c% S' S0 r/ w5 p
    tic
    ; t7 d( f0 R( n/ c7 `9 G" Ph = waitbar(0,'Please wait...');# ]  W6 c! Q. P$ r, a1 o
    for i = 0:step:360-step
    3 l; b6 D. s- {- j! K    % Delay-and-sum beamforming
    # `8 H' h* j. m5 F& b: v- C    [ DS, x1] = DelaySumURA(signal,fs,512,512,256,d,i/180*pi);
    5 h" ?3 S) o) A# ]" b    t = t+1;
      J3 Q( _4 K8 ^) i    %beamformed output energy" o: h: d1 a$ S7 a2 ~
        P(t) = DS'*DS;  T7 f, q  h# Y0 B7 c8 q1 _
        waitbar(i / length(step:360-step))7 y' ^8 F" b8 P
    end) A, ^; k2 P) u, o8 U1 j5 n4 d
    toc
    $ o& f% x8 |0 P7 J/ K/ Bclose(h)
    4 G: @7 U6 x) {" a: B; n7 L# ?5 d[m,index] = max(P);
    % }/ V8 r: Z4 U! {  qfigure,plot(0:step:360-step,P/max(P))
    ; X5 `7 l( }; G' D4 v% Qang = (index)*step; I9 n2 _" W" G0 {4 `
    # k6 u( c- h6 S: P
    程序中用的是圆阵,可以进行二维方向角扫描,不过这里为了简便就固定了俯仰角,只扫描方位角,结果如下
    7 ~, R. E# K8 V# }4 _- t4 y& L5 a5 l6 @2 `

    结果与预期相同

    PHAT加权

      与GCC-PHAT方法相同,这里也可以对幅度做归一化,只保留相位信息,使得到的峰值更明显,提高在噪声及混响环境下的性能 " ?! a) ^1 Q7 Z
      上面代码中加上这一句

    ) r# l* E2 y+ A
    %x_fft = bsxfun(@rdivide, x_fft,abs(d(1:N/2+1,));
    - K; f* A/ h0 C4 c, u测试同样的文件,结果如下
    ) h! \4 C! f1 }5 I; Z
    " m$ S* }6 s# {. I1 }0 B) f% S" _对比可以看到,PHAT加权的方法性能更好- Q" g  {- W. u% m6 }
    参考; W0 e6 \7 u# F8 C7 R, r
    1.《SRP-PHAT-A High-Accuracy, Low-Latency Technique for Talker Localization in Reverberant Environments Using Microphone Arrays》
    1 ~* a# e2 @! [& q2. 《传感器阵列波束优化设计与应用》2 C( y9 i9 F" q* }8 n2 l7 d1 \
    ————————————————" ?4 b" I" \( M2 p8 o
    版权声明:本文为CSDN博主「373955482」的原创文章,遵循CC 4.0 BY-SA版权协议,转载请附上原文出处链接及本声明。
    4 ~% _. F4 I3 d7 c6 V) t( @, V原文链接:https://blog.csdn.net/u010592995/article/details/81586504
    - O. ~) s! T/ o$ g% P$ q1 c! v+ [/ D* H; d5 M" a0 G8 S/ s
    - e0 ?' r' T4 h" G

    5 k+ h3 O6 i6 F, p4 d  i
    , z4 m5 Y- |" ~; `' g, M, B; L6 y# A# \8 m+ n0 R0 F* b2 M
    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 13:36 , Processed in 0.443054 second(s), 50 queries .

    回顶部