QQ登录

只需要一步,快速开始

 注册地址  找回密码
查看: 4612|回复: 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" F- b, i1 ~9 x4 S7 o" g; b2 E
      声源定位方法一般可分为三类,一种是基于TDOA的两步算法(two-stage algorithm),一种是基于空间谱估计如MUSIC等,还有就是基于beamforming的方法,也就是这里要介绍的可控波束响应(steered-response power),
    , J! x0 j" R0 x# f1 f5 V/ I/ g8 v: Y
    steered-response power8 c  v! u2 S+ m1 O5 U. j
      可控波束响应是利用波束形成(beamforming)的方法,对空间不同方向的声音进行增强,得到声音信号最强的方向就被认为是声源的方向。
    ; E: Z! L: A: `' S  上一篇中简单介绍了麦克风阵列的背景知识,最简单的SRP就是利用延时-累加(delay-and-sum)的方法,寻找输出能量最大的方向。
    8 ]! J. Q! K) L1 s, l% X4 |2 X  其中,语音信号为宽带信号,因此需要做宽带波束形成,这里我们在频域实现
    1 D0 J" t  ^9 _9 M# f# T+ c2 T: I5 t" A
    " K$ u, x$ @5 c0 d  B频域宽带波束形成' Z8 R4 d# Y# K( C
      频域宽带波束形成可以归类为DFT波束形成器,结构如下图
    0 H2 x. Z, t& p$ o( J; ?, Y% c- L9 c

    ! z- `5 K0 c5 q" m2 k3 W

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

    频域宽带延时累加波束形成的基本过程就是信号分帧加窗->DFT->各频点相位补偿->IDFT
    ! `. e( L; w) K/ [3 Z9 ~+ q代码实现如下

      q1 p: g0 F! E9 q+ O, V
    nction [ DS, x1] = DelaySumURA( x,fs,N,frameLength,inc,r,angle)
    & T( j- u' D5 F0 r5 B. X%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%( f! I0 D% p( j" a: s
    %frequency-domain delay-sum beamformer using circular array
    7 F' G7 H4 X& }6 P%   % T" F# W6 ^8 j1 U, ?% T; n
    %      input :, E% r: R% k: R: i+ }4 n! |3 I9 |
    %          x : input signal ,samples * channel3 B* ]  K1 K& R- W/ F3 s# U* \' D7 ^
    %          fs: sample rate1 i* C* K; ^1 \4 j& k
    %          N : fft length,frequency bin number$ Q& @$ Z4 N! n; ^8 \. y
    %frameLength : frame length,usually same as N; [# q/ L; V# S7 S
    %        inc : step increment! O1 u# f( ^2 n) u3 A
    %          r : array element radius
    9 X! {6 v4 W7 H& f7 x8 n8 V%      angle : incident angle
    0 y! m4 @7 B7 d; x2 i7 S: \5 j%6 ~( W2 M2 V1 z$ A4 D
    %     output :! n! \4 x5 e. p' j5 o0 s
    %         DS : delay-sum output
    9 i. I1 S! Z, Z- U0 O%         x1 : presteered signal,same size as x  j3 l6 }, ?3 _7 ~' w  ~
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%# Z# a# D- ]; k) P3 ~# L+ M! t

    ' M% q) j  }, G% G3 Pc = 340;
    ' e* z& o0 b0 V6 Z, tNele = size(x,2);% D' q5 O1 }) ^  p9 F1 I
    omega = zeros(frameLength,1);
    : ^) r0 z3 D5 O2 ~: Y3 EH = ones(N/2+1,Nele);
    ! H4 {1 j4 S2 A
    ' _4 j  K0 u* E/ j6 T: ytheta = 90*pi/180; %固定一个俯仰角
    , X5 c4 A+ \) p& @! H1 Hgamma = [30 90 150 210 270 330]*pi/180;%麦克风位置
    3 @: l6 n4 g9 r" d  Stao = r*sin(theta)*cos(angle(1)-gamma)/c;     %方位角 0 < angle <360, A/ k' L9 V* S2 E
    yds = zeros(length(x(:,1)),1);
    ) j( c! R  w/ ]1 wx1 = zeros(size(x));
    ' j: J! N/ l4 d1 a2 c" i, M: g  k' `; w# a% T' C* T
    % frequency bin weights, S  o" x# n% G+ p5 {
    % for k = 2:1:N/2+1, W6 d$ L2 z- k) u$ W& o- u& e
    for k = 1:1:5000*N/fs
    # |9 \+ Z7 a6 w$ R" h    omega(k) = 2*pi*(k-1)*fs/N;   
    / c3 I$ P8 G- f+ c/ k8 p& h0 h1 ]    % steering vector
    & e6 f  O1 ], Z! R5 Z  ^    H(k, = exp(-1j*omega(k)*tao);9 j; k9 Y; D% E  n$ ]4 [" E) N
    end* F$ {, Y- C$ E6 b
    * a+ ^2 z) z: Z# p$ N% H8 c
    for i = 1:inc:length(x(:,1))-frameLength  |2 C3 |' m1 w- M
    , A  D- G+ ^% @
        d = fft(bsxfun(@times, x(i:i+frameLength-1,,hamming(frameLength)));
    ! _4 S+ K2 C0 l$ F, D4 {3 w2 S' X, f
        x_fft=bsxfun(@times, d(1:N/2+1,,H);' T- }$ z1 T: m+ K" b: d0 j) }7 ]" W
    : \- U9 c$ ]% e& y; n2 F& [
        % phase transformed) o7 B9 M9 A% v# E) g# ?
        %x_fft = bsxfun(@rdivide, x_fft,abs(d(1:N/2+1,));
    8 N: }- a$ G( y6 }    yf = sum(x_fft,2);  F2 o5 B, ]; L% l# T
        Cf = [yf;conj(flipud(yf(2:N/2)))];
    , Y* D2 ~! }  e' Y9 m$ |5 p, u; D+ t0 S, d5 U6 G# u2 y
        % 恢复延时累加的信号
    5 Y1 L! _$ T5 Z1 O( ^! r    yds(i:i+frameLength-1) = yds(i:i+frameLength-1)+(ifft(Cf));
    5 t5 K8 |! x* w) @# h" E
    # t. b8 b0 S) J2 F1 S( F    % 恢复各路对齐后的信号( o9 K0 W# @* |$ @7 v
        xf  = [x_fft;conj(flipud(x_fft(2:N/2,))];$ A2 j6 y* N4 h
        x1(i:i+frameLength-1, = x1(i:i+frameLength-1,+(ifft(xf));
    . v0 a! r) B. n9 O, p6 w5 m9 k3 Jend
    - l. i' x. D  b. G$ b6 ODS = yds/Nele;  5 |8 c, O% |3 A+ T- R) s

    4 j9 l% L  H! uend- J6 S- Z, O" |0 i6 k8 A
    然后遍历各个角度重复调用这个函数,测试实际录音数据,代码如下
    ' t7 k. o8 q0 l; z5 V" Y. `& H3 e2 b2 _
    , u+ f+ `( Y) {( @) r%% SRP Estimate of Direction of Arrival at Microphone Array
    + W. g) T0 f% d6 ~% Frequency-domain delay-and-sum test
    . w1 G1 X7 D8 x/ A%  
    & \$ A! t; J# W2 G0 T* h%%
    & X4 c. L: Q1 @( k8 J
    4 D# a7 A& r6 |6 c% L2 A% x = filter(Num,1,x0);2 C$ Y# f+ ?- E" U2 m
    c = 340.0;
    ; F# G: i( [+ P- N3 _$ U1 b6 R
    8 w& w2 {1 i% i8 I9 b8 H% XMOS circular microphone array radius. X$ G3 S8 r/ Y4 T9 W
    d = 0.0420;6 q  X8 \% V  V
    % more test audio file in ../../TestAudio/ folder) E! |2 m4 J! j" d& ^
    path = '../../TestAudio/XMOS/room_mic5-2/';
    ' ]( \# J/ b: N3 T  f6 E' [5 M[s1,fs] = audioread([path,'音轨-2.wav']);( `7 Y1 Y& z, q1 y; k
    s2 = audioread([path,'音轨-3.wav']);3 D: O. m' A7 C% f
    s3 = audioread([path,'音轨-4.wav']);# h1 W2 L# S1 Q( p/ N
    s4 = audioread([path,'音轨-5.wav']);
    0 }' D. [: O1 A* W# qs5 = audioread([path,'音轨-6.wav']);* e8 p% u5 _0 w/ Y: V/ s+ D
    s6 = audioread([path,'音轨-7.wav']);3 G( u- [/ v9 v2 V
    signal = [s1,s2,s3,s4,s5,s6];
      j$ d! V7 [; F, K! x- LM = size(signal,2);" k6 ]8 {3 D% F2 F
    %%" _6 ^5 I. \2 B+ z
    t = 0;% c4 d6 h. h7 t) H. C

    . l! M4 Y8 q4 F2 i, Y% minimal searching grid, K9 ~) a$ ~$ \+ }; O: v" B+ m6 A
    step = 1;
    9 T, d+ {  @! F+ i- V" f' B/ X' M% o
    P = zeros(1,length(0:step:360-step));" Z9 B3 f% t' ]7 g+ t
    tic
    # [! z9 n& g( K# _; a0 i  ^" lh = waitbar(0,'Please wait...');
      U5 q. V, O5 ]1 q9 ifor i = 0:step:360-step: Q7 ~1 M9 W% o) x* F) n
        % Delay-and-sum beamforming; e! D6 B6 |' P! _) N6 l% t, l
        [ DS, x1] = DelaySumURA(signal,fs,512,512,256,d,i/180*pi);
    5 S  S. G% w! @7 X7 H7 d9 }- l) K    t = t+1;
    7 x" x! G9 G$ m( C    %beamformed output energy
    3 N" Q5 d& e% X9 c  @    P(t) = DS'*DS;
    / R7 U  i/ ~+ |: y% ^( o5 q& Q    waitbar(i / length(step:360-step))
      s; x& V/ L0 L# w. nend! s: c7 c: ?3 v8 o- f1 c1 r
    toc8 I) I6 j* d9 l6 L: l4 I% @) k7 y
    close(h) 0 l& w2 d3 T* R7 q
    [m,index] = max(P);/ B+ ~6 S1 F5 h: {
    figure,plot(0:step:360-step,P/max(P)): C9 `% i- l+ A7 L
    ang = (index)*step
    0 a1 c' `- G2 [2 O0 @! E& w! z" d7 A- R% t8 s8 _8 i" |
    程序中用的是圆阵,可以进行二维方向角扫描,不过这里为了简便就固定了俯仰角,只扫描方位角,结果如下
    $ z6 X9 \2 a& X% I
    ) I/ ^1 O) A$ M. n) H/ L% Q

    结果与预期相同

    PHAT加权

      与GCC-PHAT方法相同,这里也可以对幅度做归一化,只保留相位信息,使得到的峰值更明显,提高在噪声及混响环境下的性能
    4 a; x4 w5 u, _$ U4 v6 I: @  上面代码中加上这一句


    . ?8 {+ N* u( `* P' W%x_fft = bsxfun(@rdivide, x_fft,abs(d(1:N/2+1,));2 J) o# a7 l4 ~" _& j8 \
    测试同样的文件,结果如下
    4 v% `* b$ A5 a2 L0 |
      @+ q+ {1 k& n1 Z+ E, f6 B- X) Q; y+ \对比可以看到,PHAT加权的方法性能更好5 J0 }& D+ j% ^
    参考
    8 F0 w) [/ o; @2 A5 g2 C! k7 E1.《SRP-PHAT-A High-Accuracy, Low-Latency Technique for Talker Localization in Reverberant Environments Using Microphone Arrays》 $ Y% x; `' j# `/ E0 i+ `
    2. 《传感器阵列波束优化设计与应用》
    , k- d7 @9 Q- P, j# g————————————————" C* e" g% m: o: o7 i$ o7 B
    版权声明:本文为CSDN博主「373955482」的原创文章,遵循CC 4.0 BY-SA版权协议,转载请附上原文出处链接及本声明。
    8 _/ s3 y2 [- X; D& M3 E原文链接:https://blog.csdn.net/u010592995/article/details/815865048 J$ J$ x% q6 ^
    2 V4 d( k! Y4 f8 T$ ?
    ( G: h- t# [+ _
    5 r4 Z1 s0 n* f$ B: C
    - m( M- l; Q$ Q3 U0 q, O' h  X

    2 n: i/ w! N/ n( d9 Z, T' n
    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 12:39 , Processed in 0.442309 second(s), 50 queries .

    回顶部