Making a sound camera

[Marco]

FFT is mostly useful as a fast way to implement convolution, and thus cross correlation. I don't know why cross correlation didn't work for you in matlab ... but it should, look at the example in their documentation. As you noticed and as I said just having a bunch of different phases for the different microphones doesn't give you much useful information, the cross correlation gives you something useful in the time domain. As I also said, when using FFT to implement cross correlation you have to zero pad the signals to twice the length (because FFT based convolution is circular).

[hamster_nz]

Can anyone give me an explanation of why that jambalaya with FFT will give me inter-mic delays? I understand that frequency domain graph will give me a sort of instant signature of the sound but I don't understand how that'll detect that one sound arrived a few ms on this mic vs that mic...

A short answer...

If you multiply the samples in one signal by the other, then integrate over a number
of samples you can get a measure of how 'alike' to signals are.

If you want to match signals at different offsets, you have to do this at each offset.
The result is a time series that shows how much the signals match at different alignments.

Testing a 16k sample signature in a set of 16k samples at each alignment will require
268,435,456 multiply-add operations. That is a lot of math. Or O(N^2) if you like big-O complexity notation.

FFT is an effective way to do this, but with less work.

You split the signature and test data into their frequency and phase components using FFT.

You multiply the components in one signature by the ones in the other. If the frequency
is only in one of them, it will be removed (as zero times anything = zero). If the frequency is in both, they will be enhanced (and the phase changed).

You are then left with the set of amplitudes that indicate which frequencies that are common to both signals, with some phase information (remember phase information = timing information)

By running the Inverse FFT on this result, where things are 'in sync' the magnitudes of the various components add up, where they are 'out of sync' on average they cancel out.

The result is a time series that shows how much the signals match at different alignments
 - exactly same numbers (except rounding!) as found by the hard method.

However, because FFT and IFFT are each an O(N*logN) complexity this is far more practical for large values of N.

So for 16k of samples, the hard way takes 268,435,456 units of work. The FFT way takes 2*16k*4.2 = 1,101,004 units of work (however the size of a unit of work may be different between the two processes, so you can't just say it will be 256 times faster).

This also brings up another point. To get good timing information from such a process, your signature must have a good mix of frequencies in it - so a sine wave will match to another sine wave of the same frequency when they are in phase, giving ambiguous results. The signature must also have structure to it - one bit of pure random noise is much like any other bit of random noise.

This is one of the reason that dolphins use 'clicks' for echolocation and radar uses 'chirps'. This technique will not allow you to accurately pinpoint the location of that annoying hum, unless you have lots of microphones to allow you to remove the ambiguity.

[daslolo]

Alright @hamster_nz, I understand that the sound needs harmonic for the FFT to do a match, fortunately, hums have harmonics, that is what bothers (something to do with the brain latching on the pattern). Thanks to your explanation I now understand that multiplying the spectrum of the various mic would give me the matching frequencies.
I don't understand how iFFT would give me offset between signals, which I think is what I need to triangulate the source.
I don't understand neither that thing you say about multiplying phase so I tried doing it in matlab. Phase edlta is 10 degrees between x and x2, the sine wave is at 10Hz so the offset between x and x2 is 1/10*10/360=0.027s (right?)
This graph shows 2, so I'm doing something wrong somewhere.

Code: [Select]

fc=10;
fs=32*fc;
t=0:1/fs:2-1/fs;
x=.5*cos(2*pi*20*t+20*pi/180) +.1*cos(2*pi*54*t+-60*pi/180) ;%time domain signal with phase shift
x2=.5*cos(2*pi*20*t+30*pi/180) +.1*cos(2*pi*54*t+-50*pi/180) ;
N=256; %FFT size
X = 1/N*fftshift(fft(x,N));%N-point complex DFT
X2 = 1/N*fftshift(fft(x2,N));%N-point complex DFT

%phase delta
X3=X.*X2;
phase3=atan2(imag(X3),real(X3))*180/pi; %phase
df=fs/N; %frequency resolution
sampleIndex = -N/2:N/2-1; %ordered index for FFT plot
f=sampleIndex*df; %x-axis index converted to ordered frequencies
plot(f,phase3.*real(X3));

[hamster_nz]

Sorry for the large swath of code, but it is everything (although using DFT not FFT)

Code: [Select]

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define POINTS   (128)
#define FFT_SIZE (POINTS/2+1)
#define PI       (3.1415926535897932384626433)

double samples1[POINTS],       samples2[POINTS], samples1flipped[POINTS];
double spectrum1_r[FFT_SIZE], spectrum1_i[FFT_SIZE];
double spectrum2_r[FFT_SIZE], spectrum2_i[FFT_SIZE];
double spectrumX_r[FFT_SIZE], spectrumX_i[FFT_SIZE];
double samplesX[POINTS];


/**********************************************************************/
void generate_signature(double *samples) {
   int i;
   for(i = 0; i < POINTS; i++)
      samples[i] = 0;

   for(i =  0; i <  5; i++) samples[i] = -1.0;
   for(i =  5; i < 10; i++) samples[i] = 1.0;
   for(i = 10; i < 15; i++) samples[i] = -1.0;
   for(i = 15; i < 20; i++) samples[i] = 1.0;
}

/**********************************************************************/
void flip_samples(double *s1, double *s2) {
   int i;
   for(i = 0; i < POINTS; i++)
     s2[i] = s1[POINTS-1-i];
}
/**********************************************************************/
void offset_sample(double *samples1, double *samples2, int offset) {
   int i;
   for(i = 0; i < POINTS; i++)
      samples2[i] = samples1[(i-offset)%POINTS];
}

/**********************************************************************/
void complex_mult(double *r1, double *i1,
                  double *r2, double *i2,
                  double *r3, double *i3) {
   int i;
   for(i = 0; i < FFT_SIZE; i++) {
      r3[i] = (r1[i] * r2[i]) - (i1[i] * i2[i]);
      i3[i] = (r1[i] * i2[i]) + (r2[i] * i1[i]);
   }
}

/**********************************************************************/
void add_noise(double *samples, double level) {
   int i;
   for(i = 0; i < POINTS; i++) {
      samples[i] += (rand()*level)/RAND_MAX;
   }
}

/**********************************************************************/
void print_max(double *samples) {
   int    max_i = 0;
   double max = samples[0];
   int i;
   for(i = 1; i < POINTS; i++) {
      if(samples[i] > max) {
        max_i = i;
        max   = samples[i];
      }
   }

   if(max_i >= POINTS/2)
      max_i -= POINTS;
   printf("Max is at %i (%10.4f)\n", max_i, max);
}

/**********************************************************************/
void dft(double *samples, double *r, double *i) {
  int b,p;
  for(b = 0; b < FFT_SIZE; b++)
  {
    r[b] = 0.0;
    i[b] = 0.0;
    for(p = 0; p <POINTS; p++)
    {
      double angle = 2*PI*b*p/POINTS;
      r[b] += samples[p] * cos(angle);
      i[b] -= samples[p] * sin(angle);
    }
  }
}

/******************************************************/
void idft(double *r, double *i, double *samples) {
  int b,p;
  double rs[FFT_SIZE];
  double is[FFT_SIZE];

  for(b = 0; b < FFT_SIZE; b++) {
     rs[b] = r[b]/(POINTS/2);
     is[b] = -i[b]/(POINTS/2);
  }
  rs[0]          = r[0]/POINTS;
  rs[FFT_SIZE-1] = r[FFT_SIZE-1]/POINTS;

  for(p = 0; p < POINTS; p++)
  {
    for(b = 0; b < FFT_SIZE; b++)
    {
      double angle = 2*PI*b*p/POINTS;
      samples[p] += (rs[b] * cos(angle)) + (is[b] * sin(angle));
    }
  }
}
/******************************************************/
int main(int argc, char *argv[]) {
   int i;
   /* Prepare the data */
   generate_signature(samples1);
   offset_sample(samples1,samples2, 50);
   add_noise(samples2, 0.8);

   /* Flip the samples to turn convolution it to correlation */
   flip_samples(samples1,samples1flipped);

   /* Take the DFT */
   dft(samples1flipped,spectrum1_r,spectrum1_i);
   dft(samples2,spectrum2_r,spectrum2_i);

   /* Multiply the spectrums */
   complex_mult(spectrum2_r, spectrum2_i,
                spectrum1_r, spectrum1_i,
                spectrumX_r, spectrumX_i);

   /* Inverse DFT to find how well they match */
   idft(spectrumX_r,spectrumX_i,samplesX);

   /* Output the results */
   for(i = 0; i < POINTS; i++) {
      printf("%4i, %10.7f, %10.7f, %10.7f\n",i, samples1[i], samples2[i], samplesX[i]);
   }
   print_max(samplesX);
   return 0;
}


In the picture below.

- The first graph is the 'signature' being looked for.

- The second is the data to look for the signature in

- The thrid graph is the output of the FFT correlation - the peak is when the how much the second signal is like the first.  Note the peak is at the offset you need to align to.

I will be the first to admit that this is just a hack, you will need to develop and test it...

[Marco]

The signature must also have structure to it - one bit of pure random noise is much like any other bit of random noise.

It's the opposite, each bit of noise is unique. A white noise source will give a really clear spike when you cross correlate two microphones receiving it.

[ebastler]

Anyway I got the FFT running, and having two cores is sweet, so I might as well use programming to make this thing.
https://github.com/laurentopia/M5-Signal-Multimeter

What's the hardware in that picture? (The little screen plus processor device to the left?) Looks neat!

[hamster_nz]

The signature must also have structure to it - one bit of pure random noise is much like any other bit of random noise.

It's the opposite, each bit of noise is unique. A white noise source will give a really clear spike when you cross correlate two microphones receiving it.

Yes, you are right.

[daslolo]

Anyway I got the FFT running, and having two cores is sweet, so I might as well use programming to make this thing.
https://github.com/laurentopia/M5-Signal-Multimeter

What's the hardware in that picture? (The little screen plus processor device to the left?) Looks neat!

M5Stack

[daslolo]

I will be the first to admit that this is just a hack, you will need to develop and test it...

Thank you so much: it does what I need. I like that you didn't bother with FFT and I never realized how slow DFT is until now
Why is the FFT size half the sample size?
In term of cpu cost, counting 5 microphones will this need 5!=120x iFFT + 5x FFT + 5x flipped FFT?
You flip X1 before doing a complex mult, is this the equivalent to a convolution? If so, could the entire thing be done with applying a NxN kernel to both signal?

If anyone wants to try it real quick on their mCPU, here is the arduino code:

Code: [Select]

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "Adafruit_GFX.h"// Hardware-specific library
#include <MCUFRIEND_kbv.h>
MCUFRIEND_kbv tft;

#define POINTS   (128)
#define FFT_SIZE (POINTS/2+1)
#define PI       (3.1415926535897932384626433)

double samples1[POINTS], samples2[POINTS], samples1flipped[POINTS];
double spectrum1_r[FFT_SIZE], spectrum1_i[FFT_SIZE];
double spectrum2_r[FFT_SIZE], spectrum2_i[FFT_SIZE];
double spectrumX_r[FFT_SIZE], spectrumX_i[FFT_SIZE];
double samplesX[POINTS];


/**********************************************************************/
void generate_signature(double *samples)
{
int i;
for (i = 0; i < POINTS; i++)
samples[i] = 0;

for (i = 0; i < 5; i++) samples[i] = -1.0;
for (i = 5; i < 10; i++) samples[i] = 1.0;
for (i = 10; i < 15; i++) samples[i] = -1.0;
for (i = 15; i < 20; i++) samples[i] = 1.0;
}

/**********************************************************************/
void flip_samples(double *s1, double *s2)
{
int i;
for (i = 0; i < POINTS; i++)
s2[i] = s1[POINTS - 1 - i];
}
/**********************************************************************/
void offset_sample(double *samples1, double *samples2, int offset)
{
int i;
for (i = 0; i < POINTS; i++)
samples2[i] = samples1[(i - offset) % POINTS];
}

/**********************************************************************/
void complex_mult(double *r1, double *i1,
double *r2, double *i2,
double *r3, double *i3)
{
int i;
for (i = 0; i < FFT_SIZE; i++)
{
r3[i] = (r1[i] * r2[i]) - (i1[i] * i2[i]);
i3[i] = (r1[i] * i2[i]) + (r2[i] * i1[i]);
}
}

/**********************************************************************/
void add_noise(double *samples, double level)
{
int i;
for (i = 0; i < POINTS; i++)
{
samples[i] += (rand()*level) / RAND_MAX;
}
}

/**********************************************************************/
void print_max(double *samples)
{
int    max_i = 0;
double max = samples[0];
int i;
for (i = 1; i < POINTS; i++)
{
if (samples[i] > max)
{
max_i = i;
max = samples[i];
}
}

if (max_i >= POINTS / 2)
max_i -= POINTS;
printf("Max is at %i (%10.4f)\n", max_i, max);
 tft.drawLine(max_i, 50, max_i,150, 0b111111111100000);
tft.setCursor(0, 170);
tft.print("Offset = ");
 tft.print(max_i);
}

/**********************************************************************/
void dft(double *samples, double *r, double *i)
{
int b, p;
for (b = 0; b < FFT_SIZE; b++)
{
r[b] = 0.0;
i[b] = 0.0;
for (p = 0; p < POINTS; p++)
{
double angle = 2 * PI*b*p / POINTS;
r[b] += samples[p] * cos(angle);
i[b] -= samples[p] * sin(angle);
}
}
}

/******************************************************/
void idft(double *r, double *i, double *samples)
{
int b, p;
double rs[FFT_SIZE];
double is[FFT_SIZE];

for (b = 0; b < FFT_SIZE; b++)
{
rs[b] = r[b] / (POINTS / 2);
is[b] = -i[b] / (POINTS / 2);
}
rs[0] = r[0] / POINTS;
rs[FFT_SIZE - 1] = r[FFT_SIZE - 1] / POINTS;

for (p = 0; p < POINTS; p++)
{
for (b = 0; b < FFT_SIZE; b++)
{
double angle = 2 * PI*b*p / POINTS;
samples[p] += (rs[b] * cos(angle)) + (is[b] * sin(angle));
}
}
}
/******************************************************/
void setup(void)
{
Serial.begin(115200);
  uint16_t ID = tft.readID();
  tft.begin(ID);
  tft.setRotation(1);
  tft.fillScreen(0x0000);

int i;
/* Prepare the data */
generate_signature(samples1);
offset_sample(samples1, samples2, 50);
add_noise(samples2, 0.8);

 double timer = micros();

/* Flip the samples to turn convolution it to correlation */
flip_samples(samples1, samples1flipped);

/* Take the DFT */
dft(samples1flipped, spectrum1_r, spectrum1_i);
dft(samples2, spectrum2_r, spectrum2_i);

/* Multiply the spectrums */
complex_mult(spectrum2_r, spectrum2_i,
spectrum1_r, spectrum1_i,
spectrumX_r, spectrumX_i);

/* Inverse DFT to find how well they match */
idft(spectrumX_r, spectrumX_i, samplesX);

tft.setCursor(0,0);
  tft.println(micros()-timer);

for (i = 0; i < POINTS; i++)
{
tft.drawLine(i, 100 - 10 * samples1[i - 1], i, 100 - 10 * samples1[i], 0b0000000000011111);
tft.drawLine(i, 100 - 10 * samples2[i - 1], i, 100 - 10 * samples2[i], 0b0000010000011111);
tft.drawLine(i, 100 - 2 * samplesX[i - 1], i, 100 - 2 * samplesX[i],   0b1111100000000000);
}
print_max(samplesX);
}

void loop(void) {}

[hamster_nz]

Thank you so much: it does what I need. I like that you didn't bother with FFT and I never realized how slow DFT is until now

You could easily make it a bit faster by not calling SIN() and COS() so much, you can make a lookup table but it will cost you in memory. That way sin() and cos() become something like:

sin = sin_table[(index*band)%table_size];
cos = sin_table[(table_size/4+index*band)%table_size];

You have to 'flip' (actually mirror in time) the needle that you are looking for in the haystack.  My code is slightly wrong, and gives you a peak with an offset that is off by one - item 0 should not move, but items[1] & items[n-1] should swap and so on.

I am not sure you will need all 120 IFFTs. You may only need to do a handful to get enough information that you can locate the source.

[daslolo]

Oh yeah, just 20x faster 
When I get the raspberry pi tomorrow I'll get to work, for real.
Anyone using bare metal things like https://ultibo.org/?