How to calculate an azimuth of a sound source in stereo signal

Question

Let say we have a stereo recording of only one source.

What is a correct algorithm to calculate its azimuth using inter channel time difference and inter channel level difference?

Answer 1

So here's a diagram:

The signal source is at $(x, y)$, left microphone is at $(-1, 0)$, right microphone is at $(+1, 0)$. $d=2$ is the distance between the microphones. $b$ is the distance from the source to the right microphone, and $a$ is the extra distance to the left microphone, which will create additional delay and further decrease the amplitude. The grey circle represents the wavefront as it hits the right microphone.

• The delay (time difference) between channels will be $T = \frac a c$, where $c$ is the speed of sound.
• Since microphones are pressure sensors, and pressure drops with inverse distance ($1/r$ law), the signal amplitude at the left microphone will be $G = \frac{b} {a+b}$ times the signal at the right microphone. (It's a level ratio, not a level difference.)

If you know $T$ and $G$ then you can find $b$: $$b=\frac{G T c}{1-G}$$

Then you know $b$ and $a+b$ and using law of cosines and Apollonius' theorem for the median of the triangle you can find $r$ and $\theta$:

$$r = \frac 1 2 \sqrt{4b^2+4ab+2a^2-d^2}$$ $$b^2=-d r \cos(\theta) + r^2 + \frac{d^2}{4}$$ $$\theta=\arccos\left( \frac{r}{d}+\frac{d}{4r}-\frac{b^2}{dr}\right)$$

Then plug in $a$ and $b$ and simplify:

$$r = \frac 1 2 \sqrt{\frac{2(Tc)^2(G^2+1)}{(G-1)^2} -d^2}$$ $$\theta=\arccos\left( \frac{r}{d}+\frac{d}{4r}-\frac{(G T c)^2}{(1-G)^2 d r}\right)$$

If you only know $T$, you can place the source along a hyperbola (shown in green) with focal points at the microphones. If you only know $G$, then you can place the source along a circle (shown in orange). The intersection of those curves is the location of the source:

I would think that accuracy becomes smeared out when the green and orange curves meet at sharp angles.

Of course if you only care about separating sounds by location, you can ignore all this and just work in (delay, ratio)-space.

Answer 2

To complete endolith's answer:

If the source is far enough (see Far Field Assumption), i.e. $a<<b$, then $G\rightarrow 1$ and you can only estimate angle, not the distance, by using $$ \tau = d/c \cos(\theta) $$

where $d$ is the microphones spacing, $c$ is the sound speed and $\tau$ is the delay.

Answer 3

A simple way to find what you're looking for is to compute the cross-correlation of the signals. You can do this in the frequency domain: $$ Z(w) = L(w)R^{*}(w)$$

L(w) and R(w) are the frequency spectrums for the "left" and "right" microphones. Take the inverse to get $z(t)$. Then, search for the maximum value in the magnitude of $z(t)$, and the time at which this maximum occurs will give you information on the delay between two relatively similar signals.

Answer 4

I will post here the algorithm I develop as I go,.

Too many things sounds nice and easy, then it doesn't go as expected.

So I would like that everyone sees what is being done and to check and recheck everything.

And, after we are done, finally an algorithm for basic acoustic localization with 2 mics will be available.

Because for now it is nowhere to be seen. A lot of works, explanations, talk about developed algorithm, but no explicit code is there that I could find.

So, in Python:

from numpy import *
import wave
import os

def normalize (x, a=12000):
    return a*(x/float(abs(amax(x))))

def uninterleave (data):
    """Given a stereo array, return separate left and right streams

    This function converts one array representing interleaved left and
    right audio streams into separate left and right arrays.  The return
    value is a list of length two.  Input array and output arrays are all
    Numpy arrays.
    This function is taken from a SWMixer module written by Nathan Whitehead.
    """
    return reshape(data, (2, -1), 'F')

def filterout (A, low=65, high=2000):
    """
    Returns positions (indices) of all values in array A which are <= low or >= high.
    """
    fo = []
    for x in xrange(len(A)):
        if abs(A[x])<=low: fo.append(x)
        if abs(A[x])>=high: fo.append(x)
    return fo

def bandpass (x, fs=44100, nfft=22050, low=65, high=2000):
    """Puts any FFT frequency bin out of low-high range to 0."""
    freqs = fft.fftfreq(nfft, 1.0/fs)
    discarded = filterout(freqs, low, high)
    for idx in xrange(len(x)):
        try: x[discarded[idx]] = 0.0
        except: break
    return x

def phaseshift (left, right, fs=44100, startwith=100, endwith=200):
    """Returns a phase shift, in number of samples, between left and right signals.
    To be used on stereo recordings."""
    # Caution: normalize() makes in place changes too,
    # so if you are planning to use the signals in raw form later, pass them in as copies.
    left  = normalize(left)
    right = normalize(right)
    # First tried approach was with cross-correlation.
    # It should work very good when left and right are not too close to pure tones.
    #n1 = argmax(correlate(left[-20:], right[-20:], mode="same", old_behavior=0))
    #n2 = argmax(correlate(right[-20:], left[-20:], mode="same", old_behavior=0))
    # Following is an adaptation of matlab code from:
    # https://www.reddit.com/r/DSP/comments/34w3d5/fft_to_measure_relative_phase_shift_of_two_signals/
    # Take last 1024 samples of left and right, and get fs/2 (good resolution) freq bins from FFT
    # Note: left and right are now half a size of the interleaved signal, so it is really fs/2, hence deviding by 4
    leftf  = fft.fft(left,  fs/4)
    rightf = fft.fft(right, fs/4)
    # Ignore bins we are definitely not interested in, using first half of FFT to avoid negatives:
    leftf  = leftf[:len(leftf)/2]
    rightf = rightf[:len(rightf)/2]
    # Supress freqs we don't want to be checked:
    leftf  = bandpass( leftf, fs/2, fs/4, startwith, endwith)
    rightf = bandpass(rightf, fs/2, fs/4, startwith, endwith)
    # Find position of a bin with maximal energy. We assume it is a leading fundamental freq,
    # so we calculate phase shift from it.
    # We average two FFTs to ensure that the freq in question is present in both
    maxbin = argmax( (leftf.real+rightf.real)/2.0 )
    print "bin=%i;" % maxbin, # Testing shows correct pickup of maxbin for a signal containing notes C3, E3, G3
                             # It picks up dead correct or very near bins of one of the notes
    # Now we get phases from all bins.
    leftf  = unwrap(angle(leftf))    
    rightf = unwrap(angle(rightf))    
    # And calculate our phase difference, finally!
    m = leftf[maxbin]-rightf[maxbin]
    return m # or abs(m)

def azimuth (dt, w=0.15, v=334):
    """Calculates an azimuth of a source using the phase shift between two channels and a spacing between microphones.
    Where azimuth line starts in a origin of a coordinate system  which is located on half a way between 2 mics and
    its X axis passes through these two, of course.
    dt --> Time delay in seconds
    w  --> Width (spacing between mics) in meters
    v  --> Velocity (Assumed speed of sound in your environment <Sorry it isn't c :D >) in meters per second
    Returns the azimuth angle in rad.
    Use rad2deg() to convert.
    """
    return arccos((1.0*dt*v)/w)

def calculate (folder):
    """This function loads PCM wave files, file by file from a folder given by same named arg.
    The function assumes file names to be constructed as:
    <some-prefix>_<deegree-at-which-the-file-was-recorded-as-integer>.wav
    To avoid any noise introduced by a start of recording, only last half second of file will be examined.
    """
    def custom (x, y):
        # Custom sorting rule function - sort a folder by recorded degrees
        x = x.split("_")[1].split(".")[0]
        y = y.split("_")[1].split(".")[0]
        return cmp(int(x), int(y))

    ldir = os.listdir(folder)
    ldir.sort(custom)
    for x in ldir:
        wf = wave.open(os.path.join(folder, x), "r")
        fs = wf.getframerate()
        nframes = fs/2 # Load only last half a second
        wf.setpos(wf.getnframes()-nframes)
        data = wf.readframes(nframes)
        data = fromstring(data, dtype=int16)
        # Extract left and right channel
        left, right = uninterleave(data)
        # Degrees at which the signal has been recorded
        print "%3sdg:" % x[x.find("_")+1:x.find(".")],
        shift = phaseshift(left, right, 44100, 130, 170)
        distance = (abs(shift)/(fs/2.0))*334
        print "WFDist=%.3f cm;" % (distance*100),
        alpha = azimuth(shift/(fs/2.0), 0.25)
        print "alpha=%.3f dg" % rad2deg(alpha)

----------------------------------------
I have a recording of C-Major chord of a very slightly pulsing celo-like synthesized instrument.
Recorded at all 360 angles. Very precisely, as the source was stationary and the mic system was revolving on specially designed platform.
The source is 40 cm away from mic system, the mics are spaced 25 cm and a speed of sound is taken to be 334 m/s (as you can see in code).
Recording was done in echoless studio, so errors from reverberations should be minimal.
This code gives me the following output:

  0dg: bin=81; WFDist=32.620 cm; alpha=nan dg
 10dg: bin=66; WFDist=3.589 cm; alpha=81.747 dg
 20dg: bin=73; WFDist=8.537 cm; alpha=70.032 dg
 30dg: bin=75; WFDist=12.160 cm; alpha=60.895 dg
 40dg: bin=83; WFDist=1.669 cm; alpha=93.829 dg
 50dg: bin=82; WFDist=27.441 cm; alpha=nan dg
 60dg: bin=75; WFDist=0.671 cm; alpha=91.538 dg
 70dg: bin=69; WFDist=0.922 cm; alpha=87.886 dg
 80dg: bin=82; WFDist=17.577 cm; alpha=45.324 dg
 90dg: bin=83; WFDist=11.118 cm; alpha=116.406 dg
100dg: bin=73; WFDist=1.022 cm; alpha=92.342 dg
110dg: bin=82; WFDist=11.370 cm; alpha=117.053 dg
120dg: bin=83; WFDist=1.501 cm; alpha=93.442 dg
130dg: bin=83; WFDist=5.881 cm; alpha=76.393 dg
140dg: bin=66; WFDist=2.062 cm; alpha=94.731 dg
150dg: bin=72; WFDist=0.030 cm; alpha=90.070 dg
160dg: bin=82; WFDist=6.528 cm; alpha=105.136 dg
170dg: bin=81; WFDist=7.132 cm; alpha=106.576 dg
180dg: bin=83; WFDist=17.091 cm; alpha=46.873 dg
190dg: bin=83; WFDist=2.753 cm; alpha=83.677 dg
200dg: bin=71; WFDist=1.206 cm; alpha=87.235 dg
210dg: bin=83; WFDist=5.956 cm; alpha=103.784 dg
220dg: bin=72; WFDist=2.758 cm; alpha=83.665 dg
230dg: bin=78; WFDist=0.887 cm; alpha=87.967 dg
240dg: bin=84; WFDist=1.931 cm; alpha=85.571 dg
250dg: bin=79; WFDist=10.439 cm; alpha=114.680 dg
260dg: bin=69; WFDist=0.957 cm; alpha=92.193 dg
270dg: bin=66; WFDist=0.123 cm; alpha=89.717 dg
280dg: bin=82; WFDist=7.935 cm; alpha=108.506 dg
290dg: bin=82; WFDist=26.221 cm; alpha=nan dg
300dg: bin=71; WFDist=12.157 cm; alpha=119.097 dg
310dg: bin=82; WFDist=5.087 cm; alpha=78.261 dg
320dg: bin=69; WFDist=12.116 cm; alpha=61.012 dg
330dg: bin=84; WFDist=7.641 cm; alpha=107.797 dg
340dg: bin=84; WFDist=13.761 cm; alpha=56.603 dg
350dg: bin=72; WFDist=16.315 cm; alpha=130.738 dg

----------------------------------------

As you can see, very inprecise and unreliable.
Why is that? What did I do wrong?

The alphas are nan where calculated delay gives greater distance than spacing between mics.
This is not possible in theory, therefore that's a first wrong thing.
Calculated delays are mostly wrong, therefore alphas differ from the real world.
So phaseshift() function isn't totally correct.

Do not be alarmed about bins for signals after 180 degrees being a little of from those of 0-180 dg.
It is actually only correct thing here, as the system have shaders that shifts the input frequency for sound waves coming from behind.
So this system can detect, and calculate an azimuth for sources that are behind it as well as for those in front.
If only calculations are correct.


When I get the correct azimuth, I will implement whole localization as suggested by endolith.

Note for matlab lovers: Python starts from 0, [] are indexing and slicing operator, () are for callables, tupples and grouping only.