Julia实现GCC-PHAT算法

实现如下：

网络上有大量关于广义时延估计的讲解和实现，其原理很好理解，本人不在赘述，但是在实现过程中我发现虽然大家截图一样，但是代码还是有些差别；本人一一尝试后没有发现鲁棒性很好代码，很多在普通的GCC算法能够正常工作（暂时不考虑噪声），但是使用PHAT加权后计算时延总算是0，但是在论文和其他文章提到PHAT算法应该是一个正交化（白化）过程，效果应该更改好。

在经过一下午查阅资料后，我发现并不是算法本身的问题，而是编写的问题。julia实现如下，首先获取声源：

using DSP, FFTW, Makie
import GLMakie##生成chirp信号
T = 2^14
e = randn(T)
fs = 44100
c = 340
d = 0.5
t = range(0,step=1/fs, length=T)
chirp_f = LinRange(2000, 2200, T)
yy = sin.(2pi .* chirp_f .* t )
yn = yy .+ e
add_window(yn) #加窗，自己实现一下就可以了
x = yn[1:1000]
y = yn[5:1004]

编写GCC算法如下：

l1, l2 = size(x, 1), size(y, 1)
outsize = l1 + l2 - 1
max_lag  = round(Int, d / c * fs) #d为麦克风距离，c为声速，fs采样频率
nffts = nextpow(2, outsize)
vaild_index =  (nffts/2+1-max_lag):1:(nffts/2+1+max_lag)
vaild_index = convert.(Int, vaild_index)
padX1 = convert.(ComplexF64, DSP._zeropad(x, nffts))
padX2 = convert.(ComplexF64, DSP._zeropad(y, nffts))
f = plan_fft!(padX1)
f * padX1
f * padX2
G12 = padX1 .* conj(padX2)
FFTW.ifft!(G12)
ans1 = fftshift(G12)[vaild_index]
center =  max_lag+1
delay_index = findmax(ans1)[2] - center ##等于4

按照上述代码GCC结果如下，符合结果。这个时候我们如果直接按照论文上的加权方式 $1/|G_{12}|$ 加权，那么我们需要更改为如下代码：

G12 = G12 ./ (abs.(G12) + eps(1.0))
FFTW.ifft!(G12)
ans1 = fftshift(G12)[vaild_index]
delay_index = findmax(ans1)[2] - center ##等于0

这个结果明显错误，因为我们是延迟了4个 $1/f_s$ 时间，而且在网络上找到的资料也都是这样写的，我也找不到理论上的错误，后面我在一个帖子（链接放在文末）上看到也有很多人遇到同样的事情。有个人在下面留言说他有个在R语言上工作很好的算法，我在用julia尝试后发现就可以了，于是便把该代码记录下来了。

首先我们要清楚PHAT加权的方式就是白噪声，因为我们发现加权后它在频域的模全部为1，同时相位不变，所以就可以得到加权后互谱的每一项应该是 $e^{(angle*im)}$ ,其中angle就为原来的angle，所以改进后算法为：

G12 = exp.(im .*angle.(G12))
FFTW.ifft!(G12)
ans2 = fftshift(G12)[vaild_index]
delay_index = findmax(ans2)[2] - center ##等于4

最终我们得到一个比之前更健壮的算法，对比GCC和GCC-PHAT效果如图：

可以看到，使用PHTA加权后波峰更加突出。

这里本人有个疑问，就是在查阅资料时有人提到如果x和y频率不同时，GCC和GCC-PHTA都将失效，请问这个是怎么回事呢?

ps：后续过程发现可能是因为只有某一频率分量，直接相除会导致结果错误，在后续仿真我发现如果提前向信号中注入高斯白噪声两种计算方式都能够生效，而后一种因为没有出发所以鲁棒性要好些，实际可根据自己使用场景选择具体算法。

最后参考链接如下：discrete signals - GCC-PHAT (Generalized cross correlation MATLAB) - Signal Processing Stack Exchange

Julia实现GCC-PHAT算法相关推荐

[Julia语言]使用Chudnovsky 算法快速计算圆周率 Pi (π) 值
测试用的电脑是一台10年老电脑,CPU型号:E3 1230V2,3.3GHZ,4核8线程,8GB内存. 用下面的Julia程序,计算1万位的Pi值,耗时为0.26秒. 作为比较,用Julia实现的另一 ...
声源定位之GCC-PHAT算法
声源定位之GCC-PHAT算法 1.GCC-PHAT介绍 Generalized Cross Correlation-Phase Transform,GCC-PHAT 广义互相关-相位变换 Notes ...
webrtc拥塞控制算法对比-GCC vs BBR vs PCC
1.前言现有集成在webrtc中的拥塞控制算法有三种, 分别是: 谷歌自研发的gcc, 谷歌自研发的BBR算法, 斯坦福大学提出的基于机器学习凸优化的PCC算法. 本文将探讨一下三个算法的区别和优缺 ...
FastICA算法类有哪些最新发表的毕业论文呢？
一.总体简介 FastICA算法的相关文献在2004年到2020年内共计65篇,主要集中在自动化技术.计算机技术.无线电电子学.电信技术.电工技术等领域,其中期刊论文50篇.会议论文7篇.专利文献8 ...
史上最全的机器学习资料（上）
摘要: 机器学习牵涉的编程语言十分之广,包括了MATLAB.Python.Clojure.Ruby等等.为了让开发者更加广泛.深入地了解机器学习,云栖社区组织翻译了GitHub Awesome Mac ...
Linux性能优化全景指南
原文:http://33h.co/kyq4z Part1Linux性能优化1性能优化性能指标高并发和响应快对应着性能优化的两个核心指标:吞吐和延时图片来自: www.ctq6.cn 应用负载角度:直 ...
[Machine learning] 国外程序员整理的机器学习资源大全
阅读目录本文汇编了一些机器学习领域的框架.库以及软件(按编程语言排序). 1. C++ 1.1 计算机视觉 CCV -基于C语言/提供缓存/核心的机器视觉库,新颖的机器视觉库 OpenCV-它提供C ...
Linux 性能优化全景指南
大家好我是坤哥之前一些朋友觉得奇怪,说你主要做 Java 的,公号怎么时不时地也推送一些 Linux 文章,其实不管你是哪个 xx 语言的工程师,要想进阶,Linux 性能优化是必备知识,举个例子 ...
SRP-PHAT综述
文章目录 1.SRP-PHAT介绍 2.改进的SRP-PHAT算法 2.1 基于随机搜索的空间收缩快速算法 2.2 由粗到精的空间收缩快速搜索算法 2.3 随机粒子滤波快速搜索算法 2.4 搜索空间聚 ...

Julia实现GCC-PHAT算法

Julia实现GCC-PHAT算法相关推荐

最新文章

热门文章