QQ登录

只需要一步,快速开始

 注册地址  找回密码
查看: 4580|回复: 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
    & p- h! w6 N/ H8 w/ J  声源定位方法一般可分为三类,一种是基于TDOA的两步算法(two-stage algorithm),一种是基于空间谱估计如MUSIC等,还有就是基于beamforming的方法,也就是这里要介绍的可控波束响应(steered-response power),
    0 d; `$ k% \" K3 `) S4 [% |6 l" x! _- y4 E1 X" e- [' f6 y
    steered-response power. T4 x$ P( z" Z# L% H9 u
      可控波束响应是利用波束形成(beamforming)的方法,对空间不同方向的声音进行增强,得到声音信号最强的方向就被认为是声源的方向。
    6 q) ]) a8 s0 G0 G  上一篇中简单介绍了麦克风阵列的背景知识,最简单的SRP就是利用延时-累加(delay-and-sum)的方法,寻找输出能量最大的方向。
    0 N' J+ v5 V' {1 p  其中,语音信号为宽带信号,因此需要做宽带波束形成,这里我们在频域实现! W4 X0 U3 `# i- W8 [1 a
    * k( i8 C1 B& d% o; T
    频域宽带波束形成
    % ]7 s1 V" h- ~" G0 w7 _& I1 |  频域宽带波束形成可以归类为DFT波束形成器,结构如下图
    9 ?" V2 V% C' ]& w, a/ G
      z7 e) D! K4 u8 N$ |$ L; t3 o
    & _: s3 O4 T$ I  T+ M, l7 j" l

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

    频域宽带延时累加波束形成的基本过程就是信号分帧加窗->DFT->各频点相位补偿->IDFT
    3 A# h' ~! }: m1 N代码实现如下

    + ^  ?0 O+ ]% G
    nction [ DS, x1] = DelaySumURA( x,fs,N,frameLength,inc,r,angle)
    * O7 P3 z; U8 i' D) _+ \%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    & F6 ]. |; Q* J- `6 z%frequency-domain delay-sum beamformer using circular array
    . J7 n& {5 e  b3 L" s%   
    " b7 b/ Q- U' [( `: x3 ^%      input :3 _) _% @. A3 Z' G, h
    %          x : input signal ,samples * channel8 t( N4 g7 y  \7 s6 h2 D$ a+ ?
    %          fs: sample rate) T: E, k8 i6 U6 A% b4 X
    %          N : fft length,frequency bin number( o% c1 P+ h3 i+ S2 k, t5 g" N9 t
    %frameLength : frame length,usually same as N7 d: ]9 Y* }; X3 _
    %        inc : step increment
    # |% X( q4 @2 W! U1 O" ~%          r : array element radius8 S2 N  ~# x1 s2 U, u* t
    %      angle : incident angle; k. F5 Y# d) k1 B: J, {: E
    %# N# x( T1 [. J9 V/ H
    %     output :
    , v5 Z1 n* X1 g%         DS : delay-sum output
    - N7 r7 P% p5 I5 }: n8 ]%         x1 : presteered signal,same size as x
    5 M/ _0 o7 w9 e6 x5 [. W%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    / U/ I; c* s4 S/ G2 v2 H9 H* |
    / Z, o2 j5 o- ^& l0 {0 _5 c# ?- Zc = 340;
    ( H2 Z4 w0 n' C* }Nele = size(x,2);
    $ O4 t4 ]$ x1 E; G0 ?: L4 b5 s5 xomega = zeros(frameLength,1);8 ~/ ?: ~' C- Y8 F0 k
    H = ones(N/2+1,Nele);
    2 c6 S% X: w$ P$ d9 z% {) D( E7 L- u; {+ ~- X/ |: Z5 n
    theta = 90*pi/180; %固定一个俯仰角
    0 r$ u) U, N" C0 f( Hgamma = [30 90 150 210 270 330]*pi/180;%麦克风位置
    . T# s/ n1 R8 E' \/ m. ]tao = r*sin(theta)*cos(angle(1)-gamma)/c;     %方位角 0 < angle <360
    - ^) g' {6 k3 e+ Vyds = zeros(length(x(:,1)),1);: S( c( T! g1 g( a; A
    x1 = zeros(size(x));
    : m; }# w7 Z& O8 t( b- n$ z7 h0 M) X( V
    % frequency bin weights
    - s# U3 C7 y3 G% y# j9 H1 T- j; H( J% for k = 2:1:N/2+10 G8 l6 A5 B: A/ L6 O; E( M
    for k = 1:1:5000*N/fs
    / C, m: F, ]: F) X* u  e    omega(k) = 2*pi*(k-1)*fs/N;   ( j$ \! P# _; l! q
        % steering vector7 h/ r; S) z+ Y  k( g# R0 T: i
        H(k, = exp(-1j*omega(k)*tao);0 J4 B, m8 ^" z) S5 g' c
    end
    8 o5 l. W0 G( I. d6 C# L0 p" I5 w3 }; |6 U2 r
    for i = 1:inc:length(x(:,1))-frameLength! M8 {7 ^' ~/ c( U0 y

    9 E0 ?( B# j! `& H/ Z    d = fft(bsxfun(@times, x(i:i+frameLength-1,,hamming(frameLength)));
    2 Y% u5 m& u  Z' |  ]% \. U( [3 D- I' i3 }/ f. C7 z1 k
        x_fft=bsxfun(@times, d(1:N/2+1,,H);
    % H% g- w$ w$ G$ D- y3 U( M
      p6 L& v# f/ T1 o    % phase transformed4 [  D1 D) k6 ^$ D; l9 w8 y
        %x_fft = bsxfun(@rdivide, x_fft,abs(d(1:N/2+1,));
    5 {# `/ I# \% h: T7 }; x4 g    yf = sum(x_fft,2);3 h0 p' u) R0 n$ }
        Cf = [yf;conj(flipud(yf(2:N/2)))];# }, ^; q; r( S6 ?

    - i3 \% V3 R  N- O& j, W    % 恢复延时累加的信号
    - Y( A3 @- M" `) g$ B0 I    yds(i:i+frameLength-1) = yds(i:i+frameLength-1)+(ifft(Cf));' J5 J' g. s6 H" _5 C8 A
    9 h0 d/ ^! W, _0 A  M
        % 恢复各路对齐后的信号- G- W" l3 q" h) |3 z/ Q
        xf  = [x_fft;conj(flipud(x_fft(2:N/2,))];4 ?3 @: B* j* [" G/ r( D
        x1(i:i+frameLength-1, = x1(i:i+frameLength-1,+(ifft(xf));
    % v! k! \2 z, V: E3 i3 _: P, Nend- q3 t! n3 v# A' z  T( @
    DS = yds/Nele;  
      X4 G8 p- M' s6 |: K& p, P& ^+ c2 P4 y
    end
    : W  q2 P. J. X2 v( O' d! g然后遍历各个角度重复调用这个函数,测试实际录音数据,代码如下
    + _+ s  ~! D- G6 M4 o
    - X. B( x# d/ p7 Y/ \" H%% SRP Estimate of Direction of Arrival at Microphone Array
    - ], k7 p$ K) D$ l% Frequency-domain delay-and-sum test/ t3 a; ^. u' u2 D+ o* n
    %  
    ! |# P7 U# i0 q$ Q, e" ]- |, M%%5 U' U) }1 g' B) u1 ~+ b8 I
    ; r+ l' ~7 X9 f! I4 N
    % x = filter(Num,1,x0);5 K' K' C8 @( J: _
    c = 340.0;
    , ~% W. c0 A8 @* A% U( g$ T) w0 ?; x$ u  Y9 o! y' E8 L9 I
    % XMOS circular microphone array radius0 H8 ]6 i( K7 A, l3 ~9 Q  e  |
    d = 0.0420;
    . p! ^, ~8 C9 h4 t( H: x: D% more test audio file in ../../TestAudio/ folder
    ; d* a. l4 S- c7 Upath = '../../TestAudio/XMOS/room_mic5-2/';. K% L) ^6 z- Y- k
    [s1,fs] = audioread([path,'音轨-2.wav']);
    1 m. J0 _& E1 ^+ M" z0 t8 ts2 = audioread([path,'音轨-3.wav']);
    0 E9 Q( B9 e; _% {: Es3 = audioread([path,'音轨-4.wav']);7 o, P1 p5 I0 [. ~" q8 @1 H
    s4 = audioread([path,'音轨-5.wav']);
    9 J7 Y4 |, g* U+ @3 Ts5 = audioread([path,'音轨-6.wav']);
    % w3 I8 b6 }. ?- d. M+ ds6 = audioread([path,'音轨-7.wav']);
    1 A% B, l$ j. Z# C0 gsignal = [s1,s2,s3,s4,s5,s6];
    ) _" b0 z& X  e5 n  Z/ _M = size(signal,2);
    0 R: \1 A* b2 F! h, e%%2 L, r, J7 \, ?, g
    t = 0;
    & `# A1 m/ r4 W. V9 ?. R7 G% f# u# E% r
    % minimal searching grid/ H8 M" u. q/ i- u7 E% `" d# B3 [
    step = 1;
    3 E3 }# N& n: h/ X3 E9 ]& s# u, B6 D- H9 N; q- _- v; H: J, z! P: z' f
    P = zeros(1,length(0:step:360-step));  ^. s$ P7 ~! n7 {
    tic, Z8 O4 P9 V' [+ s) b+ H
    h = waitbar(0,'Please wait...');( Z/ T9 D% \- {+ Q
    for i = 0:step:360-step
    ( j: [- V0 v/ C  ?    % Delay-and-sum beamforming  e6 F# }: W. M1 s- T8 w) Q
        [ DS, x1] = DelaySumURA(signal,fs,512,512,256,d,i/180*pi);$ H: q- u, p6 T2 V3 e0 l& }2 |
        t = t+1;
    5 H. o# w4 c7 q1 u    %beamformed output energy* N# s: j, Y' D" S( Y
        P(t) = DS'*DS;3 s7 b3 @( o6 g! w. F/ Y  Y
        waitbar(i / length(step:360-step))
    * i/ _9 Z" d  `5 yend7 h& j8 W2 J: k
    toc
    8 V2 K, X1 R/ v; I# ?close(h)
    1 h& x" r3 c6 b+ Q+ F3 V/ |- P[m,index] = max(P);
    ! i3 I. v# g. pfigure,plot(0:step:360-step,P/max(P))
    ( f; N* J, n; z) ]8 s1 ~) \( @ang = (index)*step
    # _6 p* ^' [' ^- d) k& c9 s/ P  u* i0 p2 g3 ~/ @7 j) ^: v9 A3 ~
    程序中用的是圆阵,可以进行二维方向角扫描,不过这里为了简便就固定了俯仰角,只扫描方位角,结果如下6 z( ^! R2 S  Y3 H

    ) W0 B+ s! D! X" o( Q1 \9 b

    结果与预期相同

    PHAT加权

      与GCC-PHAT方法相同,这里也可以对幅度做归一化,只保留相位信息,使得到的峰值更明显,提高在噪声及混响环境下的性能 2 v( c! {3 M' V; o+ b* \9 ~
      上面代码中加上这一句

    . B& ~. n) T& t) D# B3 }8 U2 i
    %x_fft = bsxfun(@rdivide, x_fft,abs(d(1:N/2+1,));
    % n- J) _+ h. B7 [' U, R测试同样的文件,结果如下
    # x# J. k% N/ ?2 r0 \5 G# F. }! V7 }# @' t0 b1 `! v/ r
    对比可以看到,PHAT加权的方法性能更好
    ; a5 m$ E, w# i8 f# i参考. I1 l: h& A( C1 X6 @8 [- C
    1.《SRP-PHAT-A High-Accuracy, Low-Latency Technique for Talker Localization in Reverberant Environments Using Microphone Arrays》
    ; C" ~5 ]! m6 R' L; `2. 《传感器阵列波束优化设计与应用》
    / f/ Z; ^. J# x7 |+ \* P' K8 k8 Y. W————————————————0 o2 a% a1 Q/ x( E7 e
    版权声明:本文为CSDN博主「373955482」的原创文章,遵循CC 4.0 BY-SA版权协议,转载请附上原文出处链接及本声明。) x, Y: {6 p& U2 Y- t& p
    原文链接:https://blog.csdn.net/u010592995/article/details/81586504
    $ m- a  m  ?: t  G
    # D2 U* I5 `. G) U; M3 ?8 z3 C" a) a+ n+ Y6 Q+ L" _
    ( m5 ~0 g" q8 V* t
    1 |* {0 V% D9 Q& \! l
    - e2 ]4 G7 |9 G* S2 t
    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-7-26 15:40 , Processed in 0.417203 second(s), 51 queries .

    回顶部