Skip to content

Instantly share code, notes, and snippets.

@endolith
Last active January 3, 2022 15:13
  • Star 11 You must be signed in to star a gist
  • Fork 1 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save endolith/376572 to your computer and use it in GitHub Desktop.
Distance between two microphones from delay of collinear source

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()
@mzalota
Copy link

mzalota commented Feb 25, 2014

Here are the imports that I used:

from scikits.audiolab import flacread
from scipy.signal import butter, lfilter, fftconvolve
from numpy import argmax

from parabolic import *

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

@endolith
Copy link
Author

Sorry for the crappy documentation. I'll improve it. Too bad Github doesn't notify you of comments on gists. >:(

@jestern77
Copy link

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 :)

@endolith
Copy link
Author

endolith commented May 15, 2016

@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.

@selfproject94
Copy link

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)``

@endolith
Copy link
Author

endolith commented Apr 6, 2017

@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

@selfproject94
Copy link

@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. : )

@endolith
Copy link
Author

@ervex94 Yes I recorded from a stereo microphone in other software.

@naveengerold
Copy link

naveengerold commented Sep 4, 2020

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

@endolith
Copy link
Author

endolith commented Sep 4, 2020

@naveengerold

What's the error?

It finds the distance using cross-correlation of the signals

@naveengerold
Copy link

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

@endolith
Copy link
Author

endolith commented Sep 16, 2020

@naveengerold

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment