GCC-PHAT는 지연 오프셋을 추정 할 때 항상 0에서 피크

Question

내 질문은 gcc-phat를 사용하여 두 오디오 신호 간의 지연을 추정 할 때 항상 제로 오프셋을 반환한다는 것입니다. 그러나 정상적인 교차 상관 Xcorr는 좋은 결과를 제공합니다. 이유를 모르겠습니다.

아래 두 가지 경우와 동일합니다.

GCC-PHAT (일반화 된 상호 상관 MATLAB)

GCC-PHAT (일반화 된 상호 상관)는 실제 오디오 신호에서 항상 지연 = 0에서 피크

신호의 DC 부분을 제거하려고했지만 작동하지 않습니다. 어떤 제안이라도 감사합니다!

Answer 1

이 코드 는 작동하는 것 같습니다 (아래에도 포함됨).

포함하는 첫 번째 예제 링크를 파이썬 화하면 해당 코드로 올바른 작업을 수행하는 것 같습니다.

Fs = 8000
dt = 1/Fs #0.125e-3
f1 = 100
tdelay = 0.625e-03 # try different values
t3 = np.linspace(0.0,1.0,Fs)
x3 = np.cos(2*np.pi*f1*t3);
x4 = np.cos(2*np.pi*f1*(t3-tdelay));
gcc_phat(x4,x3)

수율 :

(0.375, 배열 ([-0.02514735, -0.02196879, -0.01847165, ..., -0.0303065, -0.02794452, -0.02514735]))

(1에서 실제 지연을 뺀 값).

wav 파일을 재생 한 후 (감사합니다!), 문제는 녹음 장비가 4 개 채널 모두에 대해 너무 많은 상관 노이즈가 있다는 것입니다. 그것은 0에서 GCC-PHAT 피크를 제공합니다.

더 많은 (상관되지 않은) 노이즈를 추가하여 노이즈를 제거하려고하면 더 나은 결과를 얻을 수 있습니다. -5 샘플 대 XCORR의 -4 샘플.

네 가지 플롯은 다음과 같습니다.

1. 내가보고있는 한 채널의 섹션입니다.

2. 상관 관계가없는 가우스 잡음이 데이터에 추가 될 때 GCC-PHAT 중앙 부분.

3. 데이터에 노이즈가 추가되지 않은 경우 GCC-PHAT 중앙 부분.

4. 무소음 데이터의 XCORR 결과입니다.

이를 수행하는 코드는 아래 및 GitHub에 있습니다 . 원래 데이터를 필터링하는 것이 가능할 수도 있지만 전혀 도움이되지 않습니다.

import numpy as np
import statistics
from scipy.io import wavfile
import matplotlib.pyplot as plt
from IPython.display import Audio
from scipy.signal import kaiserord, lfilter, firwin, freqz


samplerate, data = wavfile.read("Q69905.wav",'rb')
samples = np.arange(20000,30000)

fir_filter = firwin(123,0.75)
channel_1 = lfilter(fir_filter,1, data[samples,0])
channel_2 = lfilter(fir_filter,1, data[samples,1])
channel_3 = lfilter(fir_filter,1, data[samples,2])
channel_4 = lfilter(fir_filter,1, data[samples,3])

noise_1 = np.random.normal(0,1000,len(channel_1))
noise_2 = np.random.normal(0,1000,len(channel_2))
noise_3 = np.random.normal(0,1000,len(channel_3))
noise_4 = np.random.normal(0,1000,len(channel_4))

print([statistics.mean(data[:,0]),  statistics.mean(data[:,1]), statistics.mean(data[:,2]), statistics.mean(data[:,3])])

# delay, gcc = gcc_phat(data[samples,0].astype(float)+10, data[samples,2].astype(float)+12, interp=1)
delay, gcc = gcc_phat(channel_1 + noise_1, channel_3 + noise_3, interp=1)
delay_no_noise, gcc_no_noise = gcc_phat(channel_1 , channel_3 , interp=1)

plt.figure(figsize=(20,30))
plt.subplot(4, 1, 1)
plt.plot(data[samples,0])
plt.subplot(4, 1, 2)
plt.plot(np.arange(-10,10),gcc[9990:10010],'.') # [9950:10050]
plt.subplot(4, 1, 3)
plt.plot(np.arange(-10,10),gcc_no_noise[9990:10010],'.') # [9950:10050]
plt.subplot(4, 1, 4)
lags,c, line, b = plt.xcorr(channel_1,channel_3)
plt.plot(lags,c,color='r')

print('GCC-PHAT: ' + str(delay))
print('XCORR: ' +  str(lags[np.argmax(c)]))
Audio(channel_1 + noise_1, rate=44100)

---

GitHub의 코드

"""
 Estimate time delay using GCC-PHAT 
 Copyright (c) 2017 Yihui Xiong
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
 You may obtain a copy of the License at
     http://www.apache.org/licenses/LICENSE-2.0
 Unless required by applicable law or agreed to in writing, software
 distributed under the License is distributed on an "AS IS" BASIS,
 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 See the License for the specific language governing permissions and
 limitations under the License.
"""

import numpy as np


def gcc_phat(sig, refsig, fs=1, max_tau=None, interp=16):
    '''
    This function computes the offset between the signal sig and the reference signal refsig
    using the Generalized Cross Correlation - Phase Transform (GCC-PHAT)method.
    '''
    
    # make sure the length for the FFT is larger or equal than len(sig) + len(refsig)
    n = sig.shape[0] + refsig.shape[0]

    # Generalized Cross Correlation Phase Transform
    SIG = np.fft.rfft(sig, n=n)
    REFSIG = np.fft.rfft(refsig, n=n)
    R = SIG * np.conj(REFSIG)

    cc = np.fft.irfft(R / np.abs(R), n=(interp * n))

    max_shift = int(interp * n / 2)
    if max_tau:
        max_shift = np.minimum(int(interp * fs * max_tau), max_shift)

    cc = np.concatenate((cc[-max_shift:], cc[:max_shift+1]))

    # find max cross correlation index
    shift = np.argmax(np.abs(cc)) - max_shift

    tau = shift / float(interp * fs)
    
    return tau, cc


def main():
    
    refsig = np.linspace(1, 10, 10)

    for i in range(0, 10):
        sig = np.concatenate((np.linspace(0, 0, i), refsig, np.linspace(0, 0, 10 - i)))
        offset, _ = gcc_phat(sig, refsig)
        print(offset)


if __name__ == "__main__":
    main()