Create a gist now

Instantly share code, notes, and snippets.

@endolith /delay.py
Last active Apr 10, 2017

What would you like to do?
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 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

Owner

endolith commented Mar 18, 2014

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

Owner

endolith commented May 15, 2016 edited

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

ervex94 commented Apr 5, 2017

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

Owner

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

ervex94 commented Apr 10, 2017

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

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