Python四元十字阵实现麦克风声源定向

 项目采用四元十字阵麦克风阵列，运用广义互相关算法实现精确到角度的声源定向，目标定向的精度依赖于时延估计的精度，选择恰当的时延估计算法至关重要。经过分析各时延估计算法的优缺点，确定采用广义互相关算法，该算法抗混响、抗干扰能力强，时延估计精度高，得到了广泛应用。广义互相关算法由互相关算法改进而来，算法首先计算两路信号的功率谱，然后采用频域滤波器对互功率谱函数进行处理, 再经逆傅里叶变换转化为时域的互相关函数，从而降低噪声对互相关函数的影响，使互相关函数的峰值更加尖锐。废话不多说，代码如下。代码保证可直接运行无bug，需要数据集私信我。

#!/usr/bin/env python
# -*-coding:utf-8-*-
import numpy as np
import matplotlib.pyplot as plt
def Cor_Cal(npA, npB, begin):
    COV = np.dot(npA, npB)
    return COV

def Cor_CalSeries(npA, npB):  
    if len(npA) != len(npB):
        print("Error len npA != len npB.")
        return []
    lenX = len(npA)
    calUseTimeLen = np.zeros(lenX * 2)
    y = np.zeros(lenX * 3)
    y[lenX:lenX * 2] = npA
    xmove = np.zeros(lenX * 5)
    xmove[lenX * 2:lenX * 3] = npB
    for i in np.arange(lenX * 2):
        x = xmove[2 * lenX - i:5 * lenX - i] 
        calUseTimeLen[i] = Cor_Cal(x, y, i)
    calUseTimeLen[0:len(calUseTimeLen) - 1] = calUseTimeLen[1:len(calUseTimeLen)]
    return calUseTimeLen

syspath = "D:\PycharmProjects\core"
if True:
    data_1 = []
    caiyang = open('D:\PycharmProjects\data\\ceshi.txt', 'r', encoding='utf-8')
    a = caiyang.read().split(' ')

    data_1 = [c.rstrip() for c in a]
    data1_1 = []
    for i in data_1:
        i = int(hex(int(str(i), 16)), 16)
        data1_1.append(i)
    oriData1 = data1_1[0::4]
    oriData2 = data1_1[1::4]
    oriData3 = data1_1[2::4]
    oriData4 = data1_1[3::4]



#########南北
if True:
    useData = {}
    sum1 = 0
    sum2 = 0
    sum3 = 0
    sum4 = 0
    for item in oriData1:
        sum1 += item
    for item in oriData2:
        sum2 += item
    for item in oriData3:
        sum3 += item
    for item in oriData4:
        sum4 += item
    mean1 = float(sum1 / len(oriData1))
    mean2 = float(sum2 / len(oriData1))
    mean3 = float(sum3 / len(oriData1))
    mean4 = float(sum4 / len(oriData1))

    oriData1 = [i - mean1 for i in oriData1[0:24000]]
    oriData2 = [i - mean1 for i in oriData2[0:24000]]
    oriData3 = [i - mean1 for i in oriData3[0:24000]]
    oriData4 = [i - mean1 for i in oriData4[0:24000]]

    useData1 = np.append(oriData1, np.zeros(24000))
    useData2 = np.append(oriData2, np.zeros(24000))
    useData3 = np.append(oriData3, np.zeros(24000))
    useData4 = np.append(oriData4, np.zeros(24000))

    tmp1 = np.fft.fft(useData1)
    tmp2 = np.fft.fft(useData2)
    tmp3 = np.fft.fft(useData3)
    tmp4 = np.fft.fft(useData4)
    tmp1_3 = tmp1 * (tmp3.conjugate())  # tmp0领先的时间
    out1_3 = np.fft.ifft(tmp1_3)
    xt1_3 = np.arange(len(tmp1_3))
    calUseFFT1_3 = np.append(out1_3.real[24001:48000], out1_3.real[0:24000])
    tmp2_4 = tmp2 * (tmp4.conjugate())   # tmp1领先的时间
    out2_4 = np.fft.ifft(tmp2_4)
    xt2_4 = np.arange(len(tmp2_4))
    calUseFFT2_4 = np.append(out2_4.real[24001:48000], out2_4.real[0:24000])


if True:
    calUseTime1_3 = Cor_CalSeries(oriData1, oriData3)
    calUseTime2_4 = Cor_CalSeries(oriData2, oriData4)
if True:
    if len(calUseFFT1_3) > len(calUseTime1_3):
        calUseFFT1_3 = calUseFFT1_3[0:len(calUseTime1_3)]
    else:
        calUseTime1_3 = calUseTime1_3[0:len(calUseFFT1_3)]

    if len(calUseFFT2_4) > len(calUseTime2_4):
        calUseFFT2_4 = calUseFFT2_4[0:len(calUseTime2_4)]
    else:
        calUseTime2_4 = calUseTime2_4[0:len(calUseFFT2_4)]
if True:
    delta1 = (calUseFFT1_3 - calUseTime1_3)
    delta2 = (calUseFFT2_4 - calUseTime2_4)



pythonOut1_3 = np.correlate(useData1, useData3, "same")  # useData1领先的时间
pythonOut1_3 = pythonOut1_3[1:]
delta3 = (pythonOut1_3 - calUseTime1_3)
pythonOut2_4 = np.correlate(useData2, useData4, "same")  # useData1领先的时间
pythonOut2_4 = pythonOut2_4[1:]
delta4 = (pythonOut2_4 - calUseTime2_4)


xt1_3 = np.arange(len(calUseFFT1_3))
xt1_3 = xt1_3 - len(calUseFFT1_3) / 2
xt1_3 = np.arange(len(calUseTime1_3))
xt1_3 = xt1_3 - len(calUseFFT1_3) / 2
xt2_4 = np.arange(len(calUseFFT2_4))
xt2_4 = xt2_4 - len(calUseFFT2_4) / 2
xt2_4 = np.arange(len(calUseTime2_4))
xt2_4 = xt2_4 - len(calUseFFT2_4) / 2
max_idx1_3 = np.argmax(calUseTime1_3)
max_x1_3, max_y1_3 = xt1_3[max_idx1_3], calUseTime1_3[max_idx1_3]
max_idx2_4 = np.argmax(calUseTime2_4)
max_x2_4, max_y2_4 = xt2_4[max_idx2_4], calUseTime2_4[max_idx2_4]

fs = 12000
T = 1 / 12000
print("Sampling site南北:", max_x1_3)
if max_x1_3 <0:
    print("方向在1侧")
else:
    print("方向在3侧")
print("Generalized cross correlation time delay:", max_x1_3 * T)
print("Sampling site东西:", max_x2_4)
if max_x2_4 <0:
    print("方向在2侧")
else:
    print("方向在4侧")
print("Generalized cross correlation time delay:", max_x2_4 * T)