3
$\begingroup$

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?

| improve this question | |
$\endgroup$
  • $\begingroup$ Do you know the microphone separation distance and microphone directional patterns? $\endgroup$ – hotpaw2 May 12 '16 at 1:09
  • $\begingroup$ Yes, it is 25 cm and let say that they are omnidirectional. They are not exactly, but to simplify things... I would like the algo to work even when I do not know the parameters. But for now I'll do with anything that works with +-5 degrees of error. Azimuth should be 0 when total left and 180 when total right. $\endgroup$ – Dalen May 12 '16 at 12:12
  • $\begingroup$ @Dalen If by "total left" you mean "signal only in left channel", that won't ever happen. Both channels will always receive the signal, but pressure levels will be different by distance according to 1/r law. $\endgroup$ – endolith May 12 '16 at 18:32
  • $\begingroup$ @endolith : Yes, of course. And it is 1/(r**2) law called the inverse square law. How can it be of use in azimuth calculation? It can only tell how far is a source from the listener, and that is, if you know/or can guess the initial source intensity. And even then, in closed space it all goes to nothing because of echoes and reverberations. Yes, if you know the distance from each mic, you can calculate azimuth, but it's an indirect and very unreliable method. I need azimuth from inter channel time and intensity differences. $\endgroup$ – Dalen May 12 '16 at 23:02
  • $\begingroup$ Yes, I know how it should be, greater intensity on one mic suggests an inclination to that side. Greater delay on one mic suggests greater inclination to the opposite side. Nice to say and think. But how much, what's the formula? How to apply it? $\endgroup$ – Dalen May 12 '16 at 23:08
2
$\begingroup$

So here's a diagram:

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.

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:

locuses

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.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Where is the origin of the orange circle? $\endgroup$ – Dalen May 14 '16 at 14:11
  • $\begingroup$ If you know the attenuation ratio but not the delay, then you know the source is somewhere along that orange circle. If you know the delay but not the attenuation ratio, then you know the source is somewhere along (the right side of) the green hyperbola $\endgroup$ – endolith May 14 '16 at 14:17
  • $\begingroup$ To know the attenuation, I must know what was the sources initial energy. Then either 1/r or 1/r**2 can be applied. But what exactly is that orange circle, how did you get it? $\endgroup$ – Dalen May 14 '16 at 14:29
  • $\begingroup$ No, G is the amplitude ratio between the channels. If left channel has an RMS amplitude of 0.6 and right channel has an RMS amplitude of 0.8, then G is 0.6/0.8 = 0.75. I generated it using Locus command in geogebra. $\endgroup$ – endolith May 14 '16 at 14:37
1
$\begingroup$

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.

| improve this answer | |
$\endgroup$
  • $\begingroup$ You claim that alpha is equal to acos(d/w), where d is delay multiplied by speed of sound, and w a spacing between mics? Correct? If so, then the time delay (phase shift) and a spacing between mics are sufficient to calculate an azimuth under far field assumtion. True? $\endgroup$ – Dalen May 14 '16 at 0:31
  • $\begingroup$ The problem with moving the microphones very close together is that the differences between their signals becomes zero. $\endgroup$ – endolith May 14 '16 at 1:45
  • $\begingroup$ @Dalen yes that is true. $\endgroup$ – ThP May 14 '16 at 4:31
  • $\begingroup$ But it is also True that acos() domain gets from 0 to 1, i.e. from 90 degrees to zero. So I suppose we guess which quadrant it is by using intensity/pressure/level. $\endgroup$ – Dalen May 14 '16 at 14:01
  • $\begingroup$ No it goes from -1 to 1 (0-180). $\endgroup$ – ThP May 14 '16 at 14:19
0
$\begingroup$

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.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thanks, but I know how to get a time delay. Now I am having trouble because a correlation gives me enormous delays which shouldn't be greater than approx. 0.0007485 seconds for the mic spacing of 25 cm and a speed of sound being 334 m/s. I am currently using standard convolution method and I think that fast convolution will not change anything. I think my problem is that the recorded signals are too close to pure tones, so I will have to calculate real phase shift from FFT, instead of looking for a delay. If you do know how to do that, well, that will be appreciated. $\endgroup$ – Dalen May 14 '16 at 10:44
  • $\begingroup$ If the signals are the same frequency, and you have high enough resolution from your FFT, you could try looking at the angles of the frequency spectrum at the maximum frequency. $\endgroup$ – soultrane May 14 '16 at 13:25
  • $\begingroup$ You assume that freq with max energy is the most significant fundamental. Well, that might actually work. But I was thinking to take differences between more than one bin, and then take a mean value. Let say three or four bins with high energy. What do you think? $\endgroup$ – Dalen May 14 '16 at 13:37
  • $\begingroup$ I've tried both of those, and I found the maximum frequency phase difference to work well $\endgroup$ – soultrane May 14 '16 at 13:38
  • $\begingroup$ Thanks! I hope it+ll work, because with convolution I was getting enormous delays resulting in distances bigger than 1 meter. I somehow got it under control by using shorter sections for a correlation, slightly bigger than maximum i.e. the 0.0007485s*22050fs. Maximum is the distance between mics. But still, it is unreliable. $\endgroup$ – Dalen May 14 '16 at 13:55
0
$\begingroup$

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.
| improve this answer | |
$\endgroup$
  • $\begingroup$ You should post this on some other site like gist.github.com, link to it from this question, and then ask new questions when you run into specific problems $\endgroup$ – endolith May 25 '16 at 1:15
  • $\begingroup$ also you'd probably be happier using pysoundfile instead of the wave module. it does normalization, float conversion, format conversion, and interleaving for you. signal, fs = sf.read(os.path.join(folder, x)) $\endgroup$ – endolith May 25 '16 at 1:18
  • $\begingroup$ also you don't have to implement cross-correlation yourself, scipy.signal.fftconvolve(a, b[::-1]) will do it. but I guess if you're throwing away some bins then it might make sense $\endgroup$ – endolith May 25 '16 at 1:25
  • $\begingroup$ @endolith : at the end it will not be so long, just calculation stuff. Well, thanks, for pysound, but wave is great for me. scipy also implements waveread and wavewrite as in matlab, but it doesn't support all features and I wouldn't use it even if it does. Just a habit. I am not calculating cross correlation here at all. I am looking at phase difference between FFT angles. fftconvolve is great, yes, but the output is the same as for numpy.correlate() as it does the same thing just in freq domain. The difference is speed. $\endgroup$ – Dalen May 25 '16 at 18:06
  • $\begingroup$ Well you should get more accurate results if you get delay from cross-correlation of the entire signal instead of just comparing phase of the strongest frequency, since it's a multi-tone signal. With pysoundfile 6 of your lines could be combined into one: (left, right), fs = sf.read(os.path.join(folder, x), start=-nframes) :D $\endgroup$ – endolith May 25 '16 at 18:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.