QQ登录

只需要一步,快速开始

 注册地址  找回密码
查看: 4613|回复: 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( Z' J0 Z- ~1 H  声源定位方法一般可分为三类,一种是基于TDOA的两步算法(two-stage algorithm),一种是基于空间谱估计如MUSIC等,还有就是基于beamforming的方法,也就是这里要介绍的可控波束响应(steered-response power),
    7 U/ k2 {. V" m, o: n! ?) n5 U0 m8 {6 {" ~1 u  y. k& r! n
    steered-response power
    / U* ]0 R8 N* \$ V7 _" ~& y4 b  可控波束响应是利用波束形成(beamforming)的方法,对空间不同方向的声音进行增强,得到声音信号最强的方向就被认为是声源的方向。
    ' p) ]4 [( g) u- ], y8 O/ p0 U0 a$ ?3 _  上一篇中简单介绍了麦克风阵列的背景知识,最简单的SRP就是利用延时-累加(delay-and-sum)的方法,寻找输出能量最大的方向。 1 D% h; g, {: M6 ^/ M
      其中,语音信号为宽带信号,因此需要做宽带波束形成,这里我们在频域实现) T$ x1 J4 e7 V) v( u
    + x& _$ i1 |3 e+ k
    频域宽带波束形成; K( x# s/ D$ i4 @2 @& d
      频域宽带波束形成可以归类为DFT波束形成器,结构如下图 : a" q) _/ W2 ^  Z

    7 N) H! F3 _6 `6 n2 |. f: q$ R
    8 E! h7 L" U, T# h* f" ]* n

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

    频域宽带延时累加波束形成的基本过程就是信号分帧加窗->DFT->各频点相位补偿->IDFT 5 a' V+ Y( g* L8 g  @
    代码实现如下

    2 C5 v" @* j9 ?4 D$ E- T) ~
    nction [ DS, x1] = DelaySumURA( x,fs,N,frameLength,inc,r,angle)
    % N4 V) o4 K) I6 K' p9 e8 Y%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  [, t) N' M. ~2 }# I) x+ g
    %frequency-domain delay-sum beamformer using circular array' G( j" M2 n* Z
    %   
    % }# n* ~' j0 D1 o; p) p; Z0 O%      input :! o3 l" t" j9 Y! T& j$ Q- U% I
    %          x : input signal ,samples * channel; z* L& t% f1 O6 R: r3 X& H
    %          fs: sample rate
    % P) X) S+ R3 G5 \8 D%          N : fft length,frequency bin number
    : h3 g& u7 m& s3 j! c%frameLength : frame length,usually same as N- R" {( O4 {" f% C- A2 l/ b
    %        inc : step increment
    ) {# V2 m) \! j8 T# C% E1 J$ n%          r : array element radius
    0 U. A6 X0 ^, E1 r: A$ [" o/ t%      angle : incident angle. x* @& M, c) w! A; I
    %1 W9 Z& K( V7 C% o5 l/ ~
    %     output :% ]1 [# q( r" }4 h0 O3 |
    %         DS : delay-sum output
    1 Y1 q4 [6 J& q9 ^6 a8 X5 H1 S%         x1 : presteered signal,same size as x; i3 u; o& J, g' [/ q
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    3 q- P8 G/ S, z
    3 m; z- A) |0 Oc = 340;! c$ `! z" L: N
    Nele = size(x,2);# y0 b" F+ Z- j4 `% Z: |: m/ t
    omega = zeros(frameLength,1);
    3 q8 u) i' }6 T3 L; sH = ones(N/2+1,Nele);
    : ]% L* r& o$ w0 s
    5 ^4 s9 n# z/ v! S0 F6 ltheta = 90*pi/180; %固定一个俯仰角4 D+ O* Y! q: w  H
    gamma = [30 90 150 210 270 330]*pi/180;%麦克风位置
    5 |9 u7 }7 V' p7 Btao = r*sin(theta)*cos(angle(1)-gamma)/c;     %方位角 0 < angle <360
    0 w7 H6 ^9 e* Q% d: qyds = zeros(length(x(:,1)),1);
    2 T- z6 T8 X7 L4 K5 s8 W$ y- O2 i7 \x1 = zeros(size(x));( L- `, |& m) c

    3 y0 X* X2 x0 @& C: G  ~% frequency bin weights* V+ V4 `; H+ b2 ?
    % for k = 2:1:N/2+1
    5 l+ v1 ?& o  V$ Z; L. tfor k = 1:1:5000*N/fs
    1 S' I% J9 R0 E, O9 D: w( p    omega(k) = 2*pi*(k-1)*fs/N;   + W/ R$ q! w: {- B, L
        % steering vector7 d. [  N& P2 }( e
        H(k, = exp(-1j*omega(k)*tao);) _7 N8 f4 }1 `3 T8 ?
    end, X7 s  A/ ~# D/ x$ ?

    ( {2 j  q! \6 N; [% `1 V; ?3 A& ?$ zfor i = 1:inc:length(x(:,1))-frameLength
    ; y+ N) ^, p6 l0 V: @9 G% a6 a; U/ s. j8 D! v; e" q5 @$ Y
        d = fft(bsxfun(@times, x(i:i+frameLength-1,,hamming(frameLength)));' d3 k+ _6 G  B1 |* T. C% ~
    7 U( o5 z% T  z0 R
        x_fft=bsxfun(@times, d(1:N/2+1,,H);
    $ K' ]4 Z- ~& V! K
    ! s; x4 T7 J9 w" X% Q0 W% s) M& D    % phase transformed3 q4 C% R2 Y" ^. |2 I
        %x_fft = bsxfun(@rdivide, x_fft,abs(d(1:N/2+1,));  q; B, d" M7 g$ j8 I
        yf = sum(x_fft,2);( R1 E' X/ i/ y3 a
        Cf = [yf;conj(flipud(yf(2:N/2)))];; x' F( A& L; f0 I

    + T& |% {  d' B& T8 f    % 恢复延时累加的信号
    * y% S. ?( v1 J3 O  y    yds(i:i+frameLength-1) = yds(i:i+frameLength-1)+(ifft(Cf));
    6 i) b$ K3 ~: M- `7 B1 B& s6 C6 M* P5 J/ X# z4 U1 Q
        % 恢复各路对齐后的信号
    # p) E9 c- ~, V! ~$ ?) y7 {    xf  = [x_fft;conj(flipud(x_fft(2:N/2,))];
    " W' O3 r! n/ O& ]9 J. K    x1(i:i+frameLength-1, = x1(i:i+frameLength-1,+(ifft(xf));
    2 ?& z7 z& d2 c. g) ^4 Zend5 H4 K+ q  Z2 K+ ?3 ~
    DS = yds/Nele;  
    + j4 f% f- u3 j! V
    3 X& w% Q0 ^7 b: E; P0 U, Wend5 F: f/ Q- U! X. T3 i
    然后遍历各个角度重复调用这个函数,测试实际录音数据,代码如下
    ( Q. B" r$ Q. j- a1 c2 r0 k, y: d6 K4 E& ]
    / M& F/ I. V. u" e%% SRP Estimate of Direction of Arrival at Microphone Array- K( O6 P2 b8 M$ |. Y
    % Frequency-domain delay-and-sum test* ?/ A5 l, b7 H, ]8 X% w, ?
    %  9 F$ {1 Z) U1 X1 F+ N1 x
    %%+ t* k8 ?2 O* E% u' P$ h: N5 [

    " q8 R  L, ^' J% x = filter(Num,1,x0);$ G6 j; r) z3 K9 s# u; k
    c = 340.0;
    ( R" {  [$ @, y! l' ]+ V
    : o3 {; T/ E4 G, Q: w; p% XMOS circular microphone array radius1 v* o: O8 i- U  s
    d = 0.0420;, A+ M7 ~8 u1 l/ u. r
    % more test audio file in ../../TestAudio/ folder& I5 T( D" X! s+ I
    path = '../../TestAudio/XMOS/room_mic5-2/';3 S7 t; M8 S) ?5 s
    [s1,fs] = audioread([path,'音轨-2.wav']);
    9 F* @" Y& c  g/ U( ls2 = audioread([path,'音轨-3.wav']);
    % ~$ N( B% l- d0 k( d& gs3 = audioread([path,'音轨-4.wav']);
    * f$ X, w/ T  c. Fs4 = audioread([path,'音轨-5.wav']);
    7 F9 l1 Q' z+ z( }8 r& ?$ v& g4 hs5 = audioread([path,'音轨-6.wav']);) ]/ b* K- x  r% P' ^# `+ O+ w$ \' d4 J
    s6 = audioread([path,'音轨-7.wav']);9 m$ y. `; N. S* b
    signal = [s1,s2,s3,s4,s5,s6];
    / b4 k, k: c. z7 U1 C  dM = size(signal,2);
    / o# c+ f4 f, A  J6 c8 o%%
    & g; j& L9 B/ `  N$ qt = 0;
    ) u* t9 t" O7 i, {: Z8 B( ~/ _8 D; K3 W( \, n
    % minimal searching grid
    - J3 R8 p- e7 s5 L3 Mstep = 1;2 P' g, G" E) f/ B+ K( c

    6 e# ^  w$ ]0 `9 O/ YP = zeros(1,length(0:step:360-step));
    % [, ~, l+ L# s7 ]) {/ a8 t2 n. @tic2 A$ O. E9 F- Q) ^. ~7 b
    h = waitbar(0,'Please wait...');
    2 `* m/ l0 Z" J$ R) A+ C& Ofor i = 0:step:360-step: ~/ ~8 s0 g! I7 @7 }
        % Delay-and-sum beamforming9 M" h4 }! H2 [, b9 `6 \& f
        [ DS, x1] = DelaySumURA(signal,fs,512,512,256,d,i/180*pi);1 V" H! x$ S# u4 l$ y
        t = t+1;
    . C2 \, x% L  @5 [1 H6 i    %beamformed output energy9 n, {; S+ O, ?: v& N5 m4 h$ l% [5 D
        P(t) = DS'*DS;
    " f/ ~' H: q1 Q! M4 |    waitbar(i / length(step:360-step))
    # G5 O: [" L: Z: g, H2 Iend" K; e) N* o7 @7 I9 I
    toc9 F. e) \/ Y  k) X- B8 b
    close(h) 5 V' v* p. p" \1 j
    [m,index] = max(P);3 w  Z6 ]: H5 ^$ Y, e( v: i
    figure,plot(0:step:360-step,P/max(P))
    $ _0 Z- P/ W" V. j6 v; Bang = (index)*step$ Z5 W! t8 H, H% t* v7 E: a
    5 ^$ T& S: C4 r1 H. _- g9 m
    程序中用的是圆阵,可以进行二维方向角扫描,不过这里为了简便就固定了俯仰角,只扫描方位角,结果如下1 q6 G0 h; a! q& Q4 J
    " v: B8 |! k) O2 A

    结果与预期相同

    PHAT加权

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

    . N: s' [; o  Y! }2 }( n, H6 A
    %x_fft = bsxfun(@rdivide, x_fft,abs(d(1:N/2+1,));: P  n2 ?2 X* ?* G) E# M; a
    测试同样的文件,结果如下
    5 ^: b. F9 ^6 ]; a& E& d( g. d) L0 d+ O( ]
    对比可以看到,PHAT加权的方法性能更好
    1 O6 N9 s" g9 D1 V参考, `) L9 `% E; O4 x+ W
    1.《SRP-PHAT-A High-Accuracy, Low-Latency Technique for Talker Localization in Reverberant Environments Using Microphone Arrays》 8 q- W; |) m. {7 C2 \7 q! X' Z
    2. 《传感器阵列波束优化设计与应用》8 n# s+ p! f# E
    ————————————————2 H; Q+ G2 w/ E' [
    版权声明:本文为CSDN博主「373955482」的原创文章,遵循CC 4.0 BY-SA版权协议,转载请附上原文出处链接及本声明。
    6 h# \% G5 Q5 c" X原文链接:https://blog.csdn.net/u010592995/article/details/81586504: O$ N% \* ]& ?2 a7 a0 p
    1 s5 E4 V! [! T% X5 r1 b  _

    7 X7 _/ h5 K: t2 L
    8 S; w/ i# p; U0 g5 z! y' |
    * S* n/ L7 _/ a1 {" J
    2 F# p; O. I9 Z/ V4 D( x- }
    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-9-9 14:36 , Processed in 0.438819 second(s), 51 queries .

    回顶部