Example of using cross-correlation to measure the delay between two sources, and then using it to calculate distance between 2 microphones.
""" | |
Measure the distance between two microphones based on the delay of a | |
clap from a point that is collinear with the two microphones. | |
Make sure your signal isn't clipping. | |
First test: | |
2 electrets 26 cm apart, hand clap | |
Distance: -25.72 cm | |
Sub-sample distance: -26.02 cm | |
Second test: | |
Laptop's stereo mics, 6.9 cm apart, snapping fingers to one side | |
Distance: 6.79 cm | |
Sub-sample distance: 6.85 cm | |
""" | |
from __future__ import division | |
import soundfile as sf | |
from scipy.signal import butter, lfilter, fftconvolve | |
from numpy import argmax | |
# download from https://gist.github.com/endolith/255291#file-parabolic-py | |
from parabolic import parabolic | |
c = 343 # m/s = speed of sound | |
signal, fs = sf.read('stereo recording.flac') | |
# Quick high-pass filter | |
cutoff = 500 # Hz | |
B, A = butter(1, cutoff / (fs/2), btype='high') | |
signal = lfilter(B, A, signal, axis=0) | |
# Cross-correlation of the two channels (same as convolution with one reversed) | |
corr = fftconvolve(signal[:, 0], signal[::-1, 1], mode='same') | |
# Find the offset of the peak. Zero delay would produce a peak at the midpoint | |
delay = int(len(corr)/2) - argmax(corr) | |
distance = delay / fs * c | |
print("Distance: %.2f cm" % (distance * 100)) | |
# More accurate sub-sample estimation with parabolic interpolation: | |
delay = int(len(corr)/2) - parabolic(corr, argmax(corr))[0] | |
distance = delay / fs * c | |
print("Sub-sample distance: %.2f cm" % (distance * 100)) |
""" | |
Listens for impulsive sounds. When one is heard, it displays the recording of | |
the left and right channels on the top, and their cross-correlation on the | |
bottom, finds the peak using parabolic/quadratic interpolation, plots | |
the peak as a red dot and prints the measured distance between the microphones. | |
If you disable the crest factor check, and play white noise with a phone, you | |
can see the red dot move back and forth as you vary the angle to the | |
microphones. | |
""" | |
from __future__ import division, print_function | |
import numpy as np | |
from matplotlib.mlab import rms_flat | |
from matplotlib import pyplot as plt | |
import pyaudio | |
from scipy.signal import fftconvolve, butter, lfilter | |
from PyQt4.QtGui import QApplication | |
from parabolic import parabolic | |
fs = 96000 # sampling rate | |
format = pyaudio.paFloat32 # max = 1.0 | |
channels = 2 | |
chunk = int(fs/4) | |
c = 343 # m/s = speed of sound | |
def crest_factor(signal): | |
""" | |
Crest factor of a 1D signal | |
""" | |
peak = np.amax(np.absolute(signal)) | |
rms = rms_flat(signal) | |
return peak / rms | |
def callback(in_data, frame_count, time_info, status): | |
""" | |
Called on each incoming frame to process data | |
""" | |
global result | |
global result_waiting | |
if in_data: | |
print('.', end='') | |
result = np.fromstring(in_data, dtype=np.float32) | |
result = np.reshape(result, (chunk, 2)) # stereo | |
result_waiting = True | |
else: | |
print('no input') | |
return None, pyaudio.paContinue | |
# Initialize blank plots | |
plt.figure(1) | |
plt.subplot(2, 1, 1) | |
t = np.arange(chunk)/fs | |
plt_L = plt.plot(t, np.ones(chunk), 'blue')[0] | |
plt_R = plt.plot(t, -np.ones(chunk), 'red')[0] # Red = Right | |
plt.margins(0, 0.1) | |
plt.grid(True, color='0.7', linestyle='-', which='major', axis='both') | |
plt.grid(True, color='0.9', linestyle='-', which='minor', axis='both') | |
plt.title('Recording') | |
plt.xlabel('Time [seconds]') | |
plt.ylabel('Amplitude') | |
plt.subplot(2, 1, 2) | |
lags = np.arange(chunk)-chunk/2 | |
plt_corr = plt.plot(lags, np.zeros(chunk))[0] | |
plt_peak = plt.plot(0, 0, 'ro')[0] | |
plt.margins(0, 0.1) | |
plt.grid(True, color='0.7', linestyle='-', which='major', axis='both') | |
plt.grid(True, color='0.9', linestyle='-', which='minor', axis='both') | |
plt.title('Cross-correlation %.2f' % 0) | |
plt.xlabel('Lag [samples]') | |
plt.ylabel('Amplitude') | |
plt.margins(0.1, 0.1) | |
plt.xlim(-100, 100) | |
plt.ylim(-10, 10) | |
# Generate quick high-pass filter for motor hums, etc | |
cutoff = 500 # Hz | |
B, A = butter(1, cutoff / (fs/2), btype='high') | |
# Initialize global variables used by callback function | |
result = np.zeros((chunk, channels)) | |
result_waiting = False | |
p = pyaudio.PyAudio() | |
if not p.is_format_supported(fs, 0, channels, format): | |
p.terminate() | |
raise RuntimeError('Format not supported') | |
# open stream using callback (3) | |
stream = p.open(format=format, | |
channels=channels, | |
rate=fs, | |
output=False, | |
input=True, | |
frames_per_buffer=chunk, | |
stream_callback=callback) | |
print('Press Ctrl+C or close plot window to stop') | |
try: | |
stream.start_stream() | |
try: | |
while plt.fignum_exists(1): # user has not closed plot | |
if result_waiting: | |
result_waiting = False | |
# High-pass filter | |
result = lfilter(B, A, result, axis=0) | |
sig_L = result[:, 0] | |
sig_R = result[:, 1] | |
# Only update plots on impulsive sound | |
# (Disable this for continuous tracking of continuous sources, | |
# like a phone playing white noise) | |
cf = crest_factor(sig_L) | |
if cf > 18: | |
plt_L.set_data(t, sig_L) | |
plt_R.set_data(t, sig_R) | |
corr = fftconvolve(sig_R, sig_L[::-1], mode='same') | |
# Update plots | |
plt.subplot(2, 1, 1) | |
plt.title('Recording (crest factor: {:.2f})'.format(cf)) | |
plt.subplot(2, 1, 2) | |
plt_corr.set_data(lags, corr) | |
argpeak, amppeak = parabolic(corr, np.argmax(corr)) | |
plt_peak.set_data(argpeak-chunk/2, amppeak) | |
plt.ylim(np.amin(corr)*1.1, np.amax(corr)*1.1) | |
distance = (argpeak-chunk/2) / fs * c # m | |
plt.title('Cross-correlation ' | |
'(dist: {:.2f} cm)'.format(distance * 100)) | |
plt.draw() | |
# doesn't work in Spyder without this | |
# https://code.google.com/p/spyderlib/issues/detail?id=459 | |
QApplication.processEvents() | |
except KeyboardInterrupt: | |
print('\nCtrl+C: Quitting') | |
else: | |
print('\nFigure closed: Quitting') | |
finally: | |
plt.close('all') | |
stream.stop_stream() | |
stream.close() | |
p.terminate() |
Sorry for the crappy documentation. I'll improve it. Too bad Github doesn't notify you of comments on gists. >:(
HI I think I can't input any audio as I get no graph. Is there something specific to do to get pyaudio to listen to my external soundcard? I get only dots accumulating :)
@jestern77 I don't know, I would suggest doing the pyaudio examples first and getting that to work.
Also I've had more luck with pysoundfile than scikits.audiolab, I should update this.
hi endolith, thank for this example of coding but somehow my question will goes same as jestern. I just wonder how do you make stereo recording. I know the basic code which record and play, I just have no idea how does 2 microphone can be recorded . did i need to make 2 of this statement code? im sorry for my bad english # start Recording
stream =audio.open(format=FORMAT,
channels=CHANNELS,
rate=RATE, input=True,
frames_per_buffer=CHUNK)
frames = []
for i in range(0, int(RATE / CHUNK * RECORD_SECONDS)):
data = stream.read(CHUNK)
frames.append(data)``
@ervex94 I recorded the .flac file in a regular audio recording software. You can do live recording with PyAudio, go through the examples on https://people.csail.mit.edu/hubert/pyaudio/#docs
@endolith i kind of confused. Did you using another software to capture/record sound and then perform it on real time system? please correct me. : )
@ervex94 Yes I recorded from a stereo microphone in other software.
i tested delay.py worked somewhat but throws error when add large audio file, or noise wav.
@endolith how your finding the distance? any other method for determining the distance between the mic
I got this @endolith when using clear audio samples.
Distance: -21266 cm
Sub-sample distance: -21251 cm
when using large audio file and noisy audio file got this error.
corr = fftconvolve(signal[:, 0], signal[::-1, 1], mode='full')
IndexError: too many indices for array
Distance: -21266 cm
You'll have to plot the waves and the correlation and see why it's not working
IndexError: too many indices for array
Check that the shape of your arrays matches the example
Here are the imports that I used:
Please note that you need to download endolith's parabolic.py file from https://gist.github.com/endolith/255291#file-parabolic-py
I tested delay.py on my .wav file and it seemed to work... somewhat