麦克风阵列声源定位 GCC-PHAT
麦克风阵列声源定位 GCC-PHAT
麦克风阵列声源定位(一)
利用麦克风阵列可以实现声源到达方向估计(direction-of-arrival (DOA) estimation),DOA估计的其中一种方法是计算到达不同阵元间的时间差,另外一种可以看这里,这篇主要介绍经典的GCC-PHAT方法
背景
简单说明问题背景,信号模型如下图,远场平面波,二元阵列
要计算得到θθ,其实就是要求两个阵元接收到的信号时间差,现在问题变成到达时间差估计(Time-Difference-of-Arrival Estimation),因此,基于延时估计的DOA方法,其实也可以看做是分两步进行的,第一步是估计延时,第二步是计算角度,与之相对应的基于空间谱估计的DOA方法就是一步完成的。下面就分两步进行介绍
1.延时估计
1.1.互相关函数(cross-correlation function
计算y1(k)y1(k)与y2(k)y2(k)的时间差,可以计算两个信号的互相关函数,找到使互相关函数最大的值即是这两个信号的时间差
离散信号的互相关函数
R(τ)=E[x1(m)x2(m+τ)]R(τ)=E[x1(m)x2(m+τ)]
求时间差就是找到互相关函数最大时的点
D=argmaxR(n)D=argmaxR(n)
说的那么简单,那就用代码验证下
%%
% Load the chirp signal.
load chirp;
c = 340.0;
Fs = 44100;
%%d = 0.25;
N = 2;
mic = phased.OmnidirectionalMicrophoneElement;
% array = phased.URA([N,N],[0.0724,0.0418],'Element',mic);
array = phased.ULA(N,d,'Element',mic);%%
% Simulate the incoming signal using the |WidebandCollector| System
% object(TM).
arrivalAng = 42;
collector = phased.WidebandCollector('Sensor',array,'PropagationSpeed',c,...'SampleRate',Fs,'ModulatedInput',false);
signal = collector(y,arrivalAng);x1 = signal(:,1);
x2 = signal(:,2);N =length(x2);
xc = xcorr(x1,x2,'biased');
[k,ind] = max(xc);
an = acos((ind-N)/Fs*340/d)*180/pixc12 = zeros(2*N-1,1);
m = 0;
for i = -(N-1):N-1m = m+1;for t = 1:Nif 0<(i+t)&&(i+t)<=Nxc12(m) = xc12(m) + x2(t)*x1(t+i);end end
end
xc12 = xc12/N;
以上程序中的循环就是上面的定义公式,运行程序可以看到循环部分计算的互相关与直接调用matlab的xcorr结果相同(注意matlab中互相关默认没做归一化),找到互相关函数的最大值就可以得到时间差
1.2.广义互相关(generalized cross-correlation)
理论上使用上面个介绍的CCF方法就可以得到时间差,但是实际的信号会有噪声,低信噪比会导致互相关函数的峰值不够明显,这会在找极值的时候造成误差。
为了得到具有更陡峭极值的互相关函数,一般在频域使用一个加权函数来白化输入信号,这就是经典的广义互相关方法。
由维纳-辛钦定理可知,随机信号的自相关函数和功率谱密度函数服从一对傅里叶变换的关系,即x1、x2x1、x2的互功率谱可由下式计算
P(ω)=∫+∞−∞R(τ)e−jωτdτP(ω)=∫−∞+∞R(τ)e−jωτdτ
R(τ)=∫+∞−∞P(ω)ejωτdωR(τ)=∫−∞+∞P(ω)ejωτdω
这一步是把互相关函数变换到了频域,哦对,上面说到是想白化互相关函数,那就把上面第二式添加一个系数
R~(τ)=∫+∞−∞A(ω)P(ω)ejωτdωR~(τ)=∫−∞+∞A(ω)P(ω)ejωτdω
设计不同的频域系数A(ω)A(ω)对应着不同方法,这里只介绍 PHAT(phase transform)方法,即取系数如下:
A(ω)=1|P(ω)|A(ω)=1|P(ω)|
基本思想就是求时间差只需要相位信息,舍弃不相关的幅度信息以提高健壮性,可以看到当A(ω)=1A(ω)=1的情况下就是经典互相关
P(ω)P(ω)为复数,可以表示为|P(ω)|∗e−jωp|P(ω)|∗e−jωp,去掉幅度信息后,就只剩相位信息e−jωpe−jωp了,要得到相位信息,可以用P(ω)abs(P(ω))P(ω)abs(P(ω))计算,也可以直接用matlab中的angle函数计算,即angle(P(ω))angle(P(ω)),
具体得到更陡峭的峰值的理论解释如下,详情参见《麦克风阵列信号处理》P198
几行代码验证下:
x1 = [1,2,3,7,9,8,3,7]';
x2 = [4,5,6,5,4,3,8,2]';[tau,R,lag] = gccphat(x1,x2) N = length(x1)+length(x2)-1;
NFFT = 32;
P = (fft(x1,NFFT).*conj(fft(x2,NFFT)));
A = 1./abs(P);
R_est1 = fftshift(ifft(A.*P));
range = NFFT/2+1-(N-1)/2:NFFT/2+1+(N-1)/2;
R_est1 = R_est1(range);R_est2 = fftshift(ifft(exp(1i*angle(P))));
R_est2 = R_est2(range);
可以看到,三种不同写法得到的R_est1 、R_est2 与matlab自带函数gccphat计算得到的R相等。
那上面例子中的宽带语音信号,用GCC-PHAT方法得到具有陡峭峰值互相关函数,找到互相关最大时的点,结合采样频率Fs与与麦克风间距dFs与与麦克风间距d,就可以得到方向信息。频域计算互相关参考另一篇博客
2.角度计算
上面的内容计算了两个麦克风的延时,实际中假设阵列中麦克风个数为NN,则所有麦克风间两两组合共有N(N−1)/2N(N−1)/2对,记第kk个麦克风坐标为(xk,yk,zk)(xk,yk,zk),声源单位平面波传播向量u⃗ =(u,v,w)u→=(u,v,w),如果麦克风k,jk,j之间的延时为τkjτkj,则根据向量关系有下式,其中c为声速,
c∗τkj=−(xk→−xj→)∗u⃗ c∗τkj=−(xk→−xj→)∗u→
这样看起来不够直观,那就代入坐标写成标量形式如下:
c∗τkj=u∗(xk−xj)+v∗(yk−yj)+w∗(zk−zj)c∗τkj=u∗(xk−xj)+v∗(yk−yj)+w∗(zk−zj)
当有多个麦克风时,每两个麦克风就可以得到一组上式,N个麦克风就会有N∗(N−1)/2个等式N个麦克风就会有N∗(N−1)/2个等式,声源单位传播向量u⃗ =(u,v,w)u→=(u,v,w)有三个未知数,因此最少只需要三组等式,也就是三个麦克风就可以计算出声源方向,这里就先假定N=3N=3,可以得到方程组如下:
c∗τ21=u∗(x2−x1)+v∗(y2−y1)+w∗(z2−z1)c∗τ21=u∗(x2−x1)+v∗(y2−y1)+w∗(z2−z1)
c∗τ31=u∗(x3−x1)+v∗(y3−y1)+w∗(z3−z1)c∗τ31=u∗(x3−x1)+v∗(y3−y1)+w∗(z3−z1)
c∗τ23=u∗(x2−x3)+v∗(y2−y3)+w∗(z2−z3)c∗τ23=u∗(x2−x3)+v∗(y2−y3)+w∗(z2−z3)
写成矩阵形式
求出u⃗ =(u,v,w)u→=(u,v,w)后,由正余弦关系就有了角度值了
θ=acos(1w)θ=acos(1w)
α=acos(usin(acos(1w)))α=acos(usin(acos(1w)))
当麦克风数量N>3N>3时,其实所有组合信息对于角度值的计算是有冗余的,这个时候可以求出所有组合的角度值,然后利用最小二乘求出最优解,这样可以利用到所有的麦克风的信息来提高角度估计的稳定性
麦克风阵列声源定位 GCC-PHAT相关推荐
- 【声源定位】 球面散乱数据插值方法/似然估计hybrid spherical interpolation/maximum likelihood (SI/ML) 麦克风阵列声源定位
1.软件版本 MATLAB2021a 2.本算法理论知识点 球面散乱数据插值方法/似然估计SI/ML 麦克风阵列声源定位 3.算法具体理论 这个部分的程序如下所示: 这个部分理论如下所示: 本文最后的 ...
- 音视频开发(40)---麦克风阵列声源定位 GCC-PHAT
麦克风阵列声源定位 GCC-PHAT 版权声明:本文为博主原创文章,未经博主允许不得转载. https://blog.csdn.net/u010592995/article/details/79735 ...
- matlab 声源定位csdn_麦克风阵列声源定位 GCC-PHAT(一)
麦克风阵列声源定位(一) 0 a" N0 Q" t t2 l$ t) F利用麦克风阵列可以实现声源到达方向估计(direction-of-arrival (DOA) estima ...
- 基于时延法的麦克风阵列声源定位分析
文章目录 一. 关于麦克风阵列 二. 关于声源定位 三. 基于广义互相关(GCC)计算时延 四. 基于时延差的声源定位法 1. 近场模型 2. 远场模型 五. 三维空间阵列的声源定位系统实现 1. 推 ...
- 麦克风阵列声源定位四通道麦克风数据库及TDOA双曲交汇定位算法实验
麦克风阵列声源定位四通道麦克风数据库建立 四通道麦克风数据库建立物理模型的建立,来源于文献:SLoClas: A DATABASE FOR JOINT SOUND LOCALIZATION AND C ...
- 麦克风阵列声源定位 SRP-PHAT
DOA 声源定位方法一般可分为三类,一种是基于TDOA的两步算法(two-stage algorithm),一种是基于空间谱估计如MUSIC等,还有就是基于beamforming的方法,也就是这里要介 ...
- 麦克风阵列声源定位效果测试
下列图片如果不清楚可以直接访问淘宝链接,从链接中的网盘资料进行拉取.从此链接看到的购买可以跟客服说,提我可以便宜50块钱~~~ 店铺链接:首页-智能语音开发者联盟-淘宝网 产品链接:https://i ...
- 麦克风阵列声源定位解决方案
其高科技: http://www.keygotech.com/cn/solution/ssl/array/noise-source-location-based-on-mic-array 一般来说,基 ...
- matlab双麦克风阵列声源定位,双麦克风阵列拾音算法可实现空间滤波和360全向自动声源定位降噪...
不仅如此,还发掘了它的隐藏技能: 送孩子去上辅导班,忙着做笔记往往就错过正在讲的知识点,自从让他带着这只录音笔去上课,回来温习的学习效果翻倍! 爸妈想要出国旅行又怕沟通障碍,自己常因为工作忙没法陪同. ...
最新文章
- 编译器设计-解析类型
- ajax请求锁屏功能
- 自定义Windows性能监视器
- 折线 没有显示_动画折线图,你还可以试试这个图表
- c语言tcp读写二进制文件,通过TCP/IP连接发送二进制文件
- 【整理】Laravel中Eloquent ORM 关联关系的操作
- poj 1815(最小割、割集)
- 阿里云无影云电脑千万级补贴,助力广东企业居家办公
- 微信小程序点击图片放大
- 从零学习AI安全监控项目【附详细代码】
- 概率论————思维导图(上岸必备)(数字特征)
- 3、MybatisPlus
- Unity Bounds的理解
- ecshop模板基础知识
- 如何用ChemDraw Prime 绘制任意弧线箭头
- 病毒性感冒和细菌性感冒怎样区分
- 图灵联邦:一个全能、多元生态的IT技术交流社区
- Linux命令之计算器bc
- [Tomcat]配置默认访问端口及Tomcat默认访问项目
- Verilog MIPS32 CPU(六)-- MDU