Seperation of Multiple Sources.

[jay-pee]

Hi,
at first I want to thank you for writing such a nice library. I hope I can contribute to pyroomacoustics in the next months.

I started to write my master thesis in the topic of multi source localisation. I wrote a script derived from your example to see how the different algorithm performance. The scenario is to separate two white noise sources with additive uncorrelated white noise (SNR = 0dB). I was very exited to see how FRIDA performance because in your paper it separated the best. However in the simulation FRIDA had the highest error (see plot). Maybe I didn't configure the used bins right. Do you have any suggestions how to improve the performance? Is the setup maybe wrong?

Output:

SRP
Recovered azimuth: [ 72. 109.] degrees
MUSIC
Recovered azimuth: [ 70. 110.] degrees
FRIDA
Recovered azimuth: [ 35.62516827 90.06240039] degrees
TOPS
Recovered azimuth: [ 68. 112.] degrees

Here is the script:

example.py:

import numpy as np
import matplotlib.pyplot as plt
import pyroomacoustics as pra
from scipy.signal import fftconvolve

#######################
# algorithms parameters
SNR = 0.    # signal-to-noise ratio
c = 343.    # speed of sound
fs = 16000  # sampling frequency
nfft = 256  # FFT size
azimuth_list_degree = [70, 110]
# used frequency bins, dependent on the method
freq_bins = {'SRP': np.arange(30, 129),
             'MUSIC': np.arange(30, 129),
             'FRIDA': np.arange(1, 129),
             'TOPS': np.arange(1, 129)}

###########################################
# We use a UCA-C with radius 4.2 cm and 7 microphones.
# This setup is used by the Amazon Echo
R = pra.circular_2D_array([0, 0], 6, 0., 0.042)
R = np.append(R, [[0], [0]], axis=1)

######
# location of original source
num_src = len(azimuth_list_degree)
azimuth_list = np.array(azimuth_list_degree) / 180. * np.pi

x_length = (nfft // 2 + 1) * nfft
azimuth_length = azimuth_list.shape[0]
mic_signals_part = np.zeros([7, x_length, azimuth_length])
for i, azimuth in enumerate(azimuth_list):
    # propagation filter bank
    propagation_vector = -np.array([np.cos(azimuth), np.sin(azimuth)])
    delays = np.dot(R.T, propagation_vector) / c * fs  # in fractional samples
    filter_bank = pra.fractional_delay_filter_bank(delays)

    # we use a white noise signal for the source
    x = np.random.randn(x_length) / azimuth_length

    # convolve the source signal with the fractional delay filters
    # to get the microphone input signals
    mic_signals_part[:, :, i] = np.array([fftconvolve(x, filter, mode='same') for filter in filter_bank])

# Sum all signal parts to one signal
mic_signals = np.sum(mic_signals_part, axis=2)

# now add the microphone noise
for signal in mic_signals:
    signal += np.random.randn(*signal.shape) * 10**(- SNR / 20)

################################
# compute the STFT frames needed
X = np.array([
    pra.stft(signal, nfft, nfft // 2, transform=np.fft.rfft).T
    for signal in mic_signals])

##############################################
# now we can test all the algorithms available
algo_names = ['SRP', 'MUSIC', 'FRIDA', 'TOPS']
fig = plt.figure(1, figsize=(5*2, 4*2), dpi=90)
i_sps = 1
for i, algo_name in enumerate(algo_names):
    # construct the new DOA object
    # the max_four parameter is necessary for FRIDA only
    doa = pra.doa.algorithms[algo_name](R, fs, nfft, c=c, max_four=4, num_src=num_src)

    # this call here perform localization on the frames in X
    doa.locate_sources(X, freq_bins=freq_bins[algo_name])
    plt.figure(1)
    ax = fig.add_subplot(2, 2, i+1, projection='polar')
    doa.polar_plt_dirac(azimuth_ref=azimuth_list)
    plt.title(algo_name, y=1.08)

    # doa.azimuth_recon contains the reconstructed location of the source
    print(algo_name)
    print('  Recovered azimuth:', np.sort(doa.azimuth_recon) / np.pi * 180., 'degrees')
plt.subplots_adjust(top=0.9,
                    bottom=0.075,
                    left=0.05,
                    right=0.95,
                    hspace=0.4,
                    wspace=0.4)
plt.show()

in doa.py the polar_plt_dirac() function has to be adjusted to have this subplot output:

        ....
        #fig = plt.figure(figsize=(5, 4), dpi=90)
        #ax = fig.add_subplot(111, projection='polar')
        ax = plt.gca()
        ....

Thank you for your help! :)

[fakufaku]

Hi Philip,
Awesome! I really appreciate the feedback on the package, so please keep it coming.
As to FRIDA, the algorithm is a little more finicky than SRP or MUSIC to work well. In our case, we had to restrict the frequency range to higher frequencies, typically 2 to 5 kHz, to get good performance. This might vary depending on the geometry of your array.
We have the files we used to run the simulation in the paper on github too, so you can check the configuration of the algorithm there.
https://github.com/LCAV/FRIDA/blob/master/figure_doa_separation.py#L150
Please let me know if you manage to make it work, and like-wise if it doesn't :)

[jay-pee]

Thank you I will have a look into it and will keep you updated :)

In the moment I don't fully understand FRIDA and TOPS algorithm. Have to read the linked paper more intensively...

But is there another way besides try and error to get the right configuration. In example for SRP and Music you can have a look at the steered power spectrum to see which bins may be suited (See plot). Is there a method to analyse which bins are right?

In the plot you can also see another beamformer which is the SRP method with a Capon beamformer. I will push it to the repository when I fixed a few bugs and found a faster implementation.

[fakufaku]

Nice plots! And looking forward to the Capon implementation.

I don't have a systematic way of doing it. In general, you have to match the wavelength of the signal to the distances between the microphones in your array. If the wavelength is very long compared to the inter-mic separation, then there is barely any phase difference between microphones and localization will be poor, as you can observe on the bottom parts of the graphs. On the contrary, if the wavelength is much shorter than the distance, you start to get aliasing, i.e. there is a phase ambiguity. Your target is the sweet spot between the two regions.

[fakufaku]

Hi @jay-pee , did you manage to have FRIDA work ?

[jay-pee]

Hey @fakufaku sorry for my late answer. My master thesis was more about getting the full spatial spectrum over time and classify this data and do a tracking over time, so you can separate sources over time. Basically, it as a post processing of the data you get out of different DOA estimation methods. I haven't tried much but if FRIDA is aiming more on the sweet spot between low resolution and aliasing it wasn't the right DOA estimation method for me because I wanted the broadband spectrum of DOA estimations. In September, I will start as a PHD student and I hope I can contribute more then I did as a Master student. Capon is still on my list for implementation.
I think you can close this issue if you want.

[fakufaku]

@jay-pee Thanks for the feedback. It seems indeed that this might not be a good fit for FRIDA. Hope to hear from you in the future then! All the best with the start of your PhD.