From 5744fd1bbb4ceca1862a8f63d4c437176cd62f7b Mon Sep 17 00:00:00 2001
From: Baozhu Zuo <baozhu.zuo@gmail.com>
Date: Tue, 25 Dec 2018 09:36:27 +0800
Subject: [PATCH] test hardware mic phase position

baozhu@bz:/tmp$ python wave_test.py record.wav
2
['448.8125', '446.65625', '-0.1875', '330.8125', '448.375', '450.59375', '451.0']
['448.78125', '448.5625', '-0.1875', '-0.125', '448.3125', '-0.125', '-450.78125']
['0.0', '0.0', '0.0', '-0.09375', '-0.09375', '-0.09375', '0.21875']
['0.0', '0.0', '0.0', '-0.09375', '-0.09375', '-0.09375', '0.21875']
['0.0', '0.0', '0.0', '-0.09375', '-0.09375', '-0.09375', '0.21875']
['0.0', '0.0', '0.0', '-0.09375', '-0.09375', '-0.09375', '0.21875']
['0.0', '0.0', '0.0', '-0.09375', '-0.09375', '-0.09375', '0.21875']
['0.0', '0.0', '0.0', '-0.09375', '-0.09375', '-0.09375', '0.21875']
['0.0', '0.0', '0.0', '-0.09375', '-0.09375', '-0.09375', '0.21875']
['0.0', '0.0', '0.0', '-0.09375', '-0.09375', '-0.09375', '0.21875']
['0.0', '0.0', '0.0', '-0.09375', '-0.09375', '-0.09375', '0.21875']
['0.0', '0.0', '0.0', '-0.09375', '-0.09375', '-0.09375', '0.21875']
['0.0', '0.0', '0.0', '-0.09375', '-0.09375', '-0.09375', '0.21875']
['0.0', '0.0', '0.0', '-0.09375', '-0.09375', '-0.09375', '0.21875']
['0.0', '0.0', '0.0', '-0.09375', '-0.09375', '-0.09375', '0.21875']
['0.0', '0.0', '0.0', '-0.09375', '-0.09375', '-0.0625', '0.21875']
['0.0', '0.0', '0.0', '-0.09375', '-0.09375', '-0.09375', '0.21875']
['0.0', '0.0', '0.0', '-0.09375', '-0.09375', '-0.09375', '0.21875']
['0.0', '0.0', '0.0', '-0.09375', '-0.09375', '-0.09375', '0.25']
['0.0', '0.0', '0.0', '-0.09375', '-0.09375', '-0.09375', '0.25']
['0.0', '0.0', '0.0', '-0.09375', '-0.09375', '-0.09375', '0.25']
['0.0', '0.0', '0.0', '-0.09375', '-0.09375', '-0.09375', '0.25']
['0.0', '0.0', '0.0', '-0.09375', '-0.09375', '-0.09375', '0.21875']
['0.0', '0.0', '0.0', '-0.09375', '-0.09375', '-0.09375', '0.21875']
['0.0', '0.0', '0.0', '-0.09375', '-0.09375', '-0.09375', '0.21875']
['0.0', '0.0', '0.0', '-0.09375', '-0.09375', '-0.09375', '0.21875']
['0.0', '0.0', '0.0', '-0.09375', '-0.09375', '-0.09375', '0.21875']
['0.0', '0.0', '0.0', '-0.09375', '-0.09375', '-0.09375', '-0.96875']
['0.03125', '0.0', '-0.0625', '-0.09375', '-0.0625', '-0.09375', '-23.0625']
['0.0', '0.0', '-0.0625', '-0.09375', '-0.0625', '-0.09375', '-23.125']
done
48000
---
 tools/phase_test.py | 83 +++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 83 insertions(+)
 create mode 100644 tools/phase_test.py

diff --git a/tools/phase_test.py b/tools/phase_test.py
new file mode 100644
index 0000000..a099f2b
--- /dev/null
+++ b/tools/phase_test.py
@@ -0,0 +1,83 @@
+import sys
+import wave
+import numpy as np
+
+
+if len(sys.argv) != 2:
+    print('Usage: {} multi.wav'.format(sys.argv[0]))
+    sys.exit(1)
+
+
+multi = wave.open(sys.argv[1], 'rb')
+rate = multi.getframerate()
+channels = multi.getnchannels()
+
+if channels <= 1:
+    sys.exit(1)
+
+N = rate 
+
+window = np.hanning(N)
+
+interp = 4*8
+max_offset = int(rate * 0.1 / 340 * interp)
+
+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)
+    #R /= np.abs(R)
+
+    cc = np.fft.irfft(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
+
+
+print(multi.getsampwidth())
+
+while True:
+    data = multi.readframes(N)
+
+    if len(data) != multi.getsampwidth() * N * channels:
+        print("done")
+        break
+
+    if multi.getsampwidth() == 2:
+        data = np.fromstring(data, dtype='int16')
+    else:
+        data = np.fromstring(data, dtype='int32')
+    ref_buf = data[0::channels]
+
+    offsets = []
+    for ch in range(1, channels):
+        sig_buf = data[ch::channels]
+        tau, _ = gcc_phat(sig_buf * window, ref_buf * window, fs=1, max_tau=max_offset, interp=interp)
+        # tau, _ = gcc_phat(sig_buf, ref_buf, fs=rate, max_tau=1)
+
+        offsets.append(tau)
+
+    print(offsets)
+
+print(multi.getframerate())
+
+multi.close()