Instantly share code, notes, and snippets.

# endolith/delay.py Last active Sep 27, 2019

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)) 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') plt_R = plt.plot(t, -np.ones(chunk), 'red') # 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)) plt_peak = plt.plot(0, 0, 'ro') 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 Author

### 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. >:(

### jestern77 commented Jan 15, 2015

 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 Author

### 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 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

### 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. : )
Owner Author

### endolith commented Mar 20, 2018

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