Spectrogram and FFT using OpenGL

Hi,

I have read the Sonographic Sound Processing page at http://designingsound.org/2013/04/24/sonographic-sound-processing/ and think to use OpenGL fragment shaders for to implement the conversion of .wav file to a spectrogram

The only thing that I think to modify compared to this page is to treat chunks of sound by rows instead colums, cf. handle each temporal chunk of sound on a different line instead on a different column

I think to directly map the PCM/ADPCM 8/4 bits data into the texture image with the stereo/quadriphonic channels mapped into R,G,B and A colors channels and make the PCM/ADPCM to float conversion directly into the fragment shader for to have a very memory space efficient sheme
(the windowing function has of course to be computed into the fragment shader for to be too really memory efficient)

How can I make the shader for to compute the FFT on a entire line, cf. a chunk of about 1/44 milliseconds = one line on a 1024 texture image , with interleaved Left/Right (or Left Front/Left Back/Right Front/Right Back) channels stored on a PCM format (or ADPCM for a quadriphonic sound), and/or exist an OpenGL’s shader library that can already handle this ?

My problem is not about the handling of stereo/quadriphonic, I can basically handle this with the copy of the left and right stéréo channels to the Left Front/Back and Right Front/Back quadriphonics channels or the PCM conversion than can easily to be handled using only an external paletted texture access

It is more about the ADPCM handling that have to compute the final channel value relative to the previous channel sample value and the computation of the Fourrier Transform that need to access alls values of a channel on a line for to compute each of the 1024 Fourrier Transform coefficients for this line

FFTs are a fairly common application of shaders (either compute shaders or fragment shaders with render-to-texture). But you’d normally use log[sub]2/sub passes for N samples. You wouldn’t try to compute a FFT in a single call because you need all of the values in each pass to be computed before starting the next pass.

Thank, GClements

You say “FFTs are a fairly common application of shaders”

=> where can I find somes complete examples, please ?

I have googled about this but don’t find something that is fully fonctionnal :frowning:

I have begin to test the Fast Fourrier Transform from http://paulbourke.net/miscellaneous/dft/ for to understand how the FFT wotk but I don’t understand results :frowning:

The real part is only 0.5f and the imaginary part is 0.0f at the wanted frequency when I use the simple signal(angle) = cos( angle )

The real part is 1.0f and the imaginary part is always 0.0f at the wanted frequency when I use the signal(angle) = cos(angle + i sin(angle)

What sort of “Fourrier transform” other than the DFT/FFT have I to handle if I want 1.0f for the real part and 1.0f for the imaginary part at the wanted frequency when I use signal(angle) = cos(angle) + i sin(angle) ???

Here are the content of files I use :

  1. fft.h

#ifndef _FFT_H_
#define _FFT_H_
/*
   This computes an in-place complex-to-complex DFT 
 x and y are the real and imaginary arrays of 2^m points.
   dir =  1 gives forward transform
   dir = -1 gives reverse transform 
*/
int DFT( int dir, int m, double *x1, double *y1 );

/*
   This computes an in-place complex-to-complex FFT 
   x and y are the real and imaginary arrays of 2^m points.
   dir =  1 gives forward transform
   dir = -1 gives reverse transform 
*/
int FFT( int dir, long m, double *x, double *y );


#endif /* _FFT_H_ */

fft.c


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

#define FALSE 0
#define TRUE  1

/*
   Direct fourier transform
*/
int DFT( int dir, int m, double *x1, double *y1 )
{
   long i,k;
   double arg;
   double cosarg,sinarg;
   double *x2=NULL,*y2=NULL;

   x2 = malloc( m * sizeof(double));
   y2 = malloc( m * sizeof(double));

   if (x2 == NULL || y2 == NULL)
      return(FALSE);

   for ( i = 0 ; i < m ; i++ ) {
      x2[i] = 0;
      y2[i] = 0;
      arg = - dir * 2.0 * 3.141592654 * (double)i / (double)m;
      for ( k = 0 ; k < m ; k++ ) {
         cosarg = cos(k * arg);
         sinarg = sin(k * arg);
         x2[i] += (x1[k] * cosarg - y1[k] * sinarg);
         y2[i] += (x1[k] * sinarg + y1[k] * cosarg);
      }
   }

   /* Copy the data back */
   if (dir == 1) {
      for ( i = 0 ; i < m; i++ ) {
         x1[i] = x2[i] / (double)m;
         y1[i] = y2[i] / (double)m;
      }
   } else {
      for ( i = 0 ; i < m ; i++ ) {
         x1[i] = x2[i];
         y1[i] = y2[i];
      }
   }

   free(x2);
   free(y2);
   return(TRUE);
}


/*
   This computes an in-place complex-to-complex FFT 
   x and y are the real and imaginary arrays of 2^m points.
   dir =  1 gives forward transform
   dir = -1 gives reverse transform 
*/
int FFT( int dir, int m, double *x, double *y )
{
   long n,i,i1,j,k,i2,l,l1,l2;
   double c1,c2,tx,ty,t1,t2,u1,u2,z;

   /* Calculate the number of points */
   n = 1;
   for ( i = 0 ; i < m ; i++ ) 
      n *= 2;

   // printf("n= %ld 
", n);

   /* Do the bit reversal */
   i2 = n >> 1;
   j = 0;
   for ( i = 0 ; i < n -1 ; i++ ) {
      if (i < j) {
         tx = x[i];
         ty = y[i];
         x[i] = x[j];
         y[i] = y[j];
         x[j] = tx;
         y[j] = ty;
      }
      k = i2;
      while  (k <= j ) {
         j -= k;
         k >>= 1;
      }
      j += k;
   }

   // printf("/* Compute the FFT */ 
");

   c1 = -1.0; 
   c2 = 0.0;
   l2 = 1;
   for (l=0;l<m;l++) {
      l1 = l2;
      l2 <<= 1;
      u1 = 1.0; 
      u2 = 0.0;
      for (j=0;j<l1;j++) {
         for (i=j;i<n;i+=l2) {
            i1 = i + l1;
            t1 = u1 * x[i1] - u2 * y[i1];
            t2 = u1 * y[i1] + u2 * x[i1];
            x[i1] = x[i] - t1; 
            y[i1] = y[i] - t2;
            x[i] += t1;
            y[i] += t2;
         }
         z =  u1 * c1 - u2 * c2;
         u2 = u1 * c2 + u2 * c1;
         u1 = z;
      }
      c2 = sqrt((1.0 - c1) / 2.0);
      if (dir == 1) 
         c2 = -c2;
      c1 = sqrt((1.0 + c1) / 2.0);
   }

   /* Scaling for forward transform */
   if (dir == 1) {
      for (i=0;i<n;i++) {
         x[i] /= n;
         y[i] /= n;
      }
   }
   
   return(TRUE);
}

test.c


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

#include "fft.h"

#define EPSILON 0.001f


int size = 16;
int size2 = 4;

double *reals;
double *images;

void CreateArrays( int n )
{
    reals  = malloc( (n+1) * sizeof(double) );
    images = malloc( (n+1) * sizeof(double) );
}

void FreeArrays()
{
    free(reals);
    free(images);
}

void InitArraysCos( int n, float freq )
{
    int i;
    
    double angle, delta;

    angle = 0.0f;
    delta = 2.0f * M_PI / (float)(n);

    for( i =  0 ; i < n ; i++ )
    {
        reals[i]  = cos( angle * freq );
        images[i] = 0.0f;

        angle += delta;
    }
}

void InitArraysCosSin( int n, float freq )
{
    int i;
    
    double angle, delta;

    angle = 0.0f;
    delta = 2.0f * M_PI / (float)(n);

    for( i =  0 ; i < n ; i++ )
    {
        reals[i]  = cos( angle * freq );
        images[i] = sin( angle * freq );

        angle += delta;
    }
}
 

void PrintArrays( int n, bool zeroes )
{
    int i;

    for( i = 0 ; i < n ; i++ )
    {
        if( (fabs(reals[i]) > EPSILON ) || (fabs(images[i]) > EPSILON ) || ( zeroes == true ) )
        {
            printf("    Freq[%d] = (%f , %f) 
", i, reals[i], images[i] );
        }
    }
}  

int main( int argc , char **argv)
{
    float freq;

    printf("FFT v0.1 by Cyclone 
");

    CreateArrays(size);

    for( freq = 1 ; freq < size / 2 ; freq++)
    {
        printf("

Init Arrays signal(angle) = cos(angle) for freq %1f Hz ... 

", freq);
        InitArraysCos(size, freq);
        // PrintArrays(size, true);

        printf("

Compute DFT  for freq %1f Hz ... 

", freq);
        InitArraysCos(size, freq);
        DFT(1, size, reals, images);
        PrintArrays(size / 2, false);

        printf("

Compute FFT for freq %1f Hz... 

", freq);
        InitArraysCos(size, freq);
        FFT(1, size2, reals, images);
         PrintArrays(size / 2 , false);

        printf("

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

");

        printf("

Init Arrays signal(angle) = cos(angle) + i sin(angle) for freq %1f Hz ... 

", freq);
        InitArraysCosSin(size, freq);
        // PrintArrays(size, true);

        printf("

Compute DFT  for freq %1f Hz ... 

", freq);
        InitArraysCosSin(size, freq);
        DFT(1, size, reals, images);
        PrintArrays(size / 2, false);

        printf("

Compute FFT for freq %1f Hz... 

", freq);
        InitArraysCosSin(size, freq);
        FFT(1, size2, reals, images);
         PrintArrays(size / 2 , false);

        printf("

    ======================= 

");

    }


    FreeArrays();

    return 0;
}

makefile


all : fft.o test

fft.o : fft.c
    gcc -c fft.c -o fft.o

test : test.c fft.o
    gcc test.c fft.o -o test -lm

clean :
    rm fft.o test

and the result of the execution of ./test


FFT v0.1 by Cyclone 


Init Arrays signal(angle) = cos(angle) for freq 1.000000 Hz ... 



Compute DFT  for freq 1.000000 Hz ... 

    Freq[1] = (0.500000 , -0.000000) 


Compute FFT for freq 1.000000 Hz... 

    Freq[1] = (0.500000 , -0.000000) 


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



Init Arrays signal(angle) = cos(angle) + i sin(angle) for freq 1.000000 Hz ... 



Compute DFT  for freq 1.000000 Hz ... 

    Freq[1] = (1.000000 , -0.000000) 


Compute FFT for freq 1.000000 Hz... 

    Freq[1] = (1.000000 , 0.000000) 


    ======================= 



Init Arrays signal(angle) = cos(angle) for freq 2.000000 Hz ... 



Compute DFT  for freq 2.000000 Hz ... 

    Freq[2] = (0.500000 , -0.000000) 


Compute FFT for freq 2.000000 Hz... 

    Freq[2] = (0.500000 , 0.000000) 


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



Init Arrays signal(angle) = cos(angle) + i sin(angle) for freq 2.000000 Hz ... 



Compute DFT  for freq 2.000000 Hz ... 

    Freq[2] = (1.000000 , -0.000000) 


Compute FFT for freq 2.000000 Hz... 

    Freq[2] = (1.000000 , 0.000000) 


    ======================= 



Init Arrays signal(angle) = cos(angle) for freq 3.000000 Hz ... 



Compute DFT  for freq 3.000000 Hz ... 

    Freq[3] = (0.500000 , -0.000000) 


Compute FFT for freq 3.000000 Hz... 

    Freq[3] = (0.500000 , 0.000000) 


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



Init Arrays signal(angle) = cos(angle) + i sin(angle) for freq 3.000000 Hz ... 



Compute DFT  for freq 3.000000 Hz ... 

    Freq[3] = (1.000000 , -0.000000) 


Compute FFT for freq 3.000000 Hz... 

    Freq[3] = (1.000000 , -0.000000) 


    ======================= 



Init Arrays signal(angle) = cos(angle) for freq 4.000000 Hz ... 



Compute DFT  for freq 4.000000 Hz ... 

    Freq[4] = (0.500000 , -0.000000) 


Compute FFT for freq 4.000000 Hz... 

    Freq[4] = (0.500000 , 0.000000) 


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



Init Arrays signal(angle) = cos(angle) + i sin(angle) for freq 4.000000 Hz ... 



Compute DFT  for freq 4.000000 Hz ... 

    Freq[4] = (1.000000 , -0.000000) 


Compute FFT for freq 4.000000 Hz... 

    Freq[4] = (1.000000 , 0.000000) 


    ======================= 



Init Arrays signal(angle) = cos(angle) for freq 5.000000 Hz ... 



Compute DFT  for freq 5.000000 Hz ... 

    Freq[5] = (0.500000 , -0.000000) 


Compute FFT for freq 5.000000 Hz... 

    Freq[5] = (0.500000 , 0.000000) 


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



Init Arrays signal(angle) = cos(angle) + i sin(angle) for freq 5.000000 Hz ... 



Compute DFT  for freq 5.000000 Hz ... 

    Freq[5] = (1.000000 , -0.000000) 


Compute FFT for freq 5.000000 Hz... 

    Freq[5] = (1.000000 , 0.000000) 


    ======================= 



Init Arrays signal(angle) = cos(angle) for freq 6.000000 Hz ... 



Compute DFT  for freq 6.000000 Hz ... 

    Freq[6] = (0.500000 , -0.000000) 


Compute FFT for freq 6.000000 Hz... 

    Freq[6] = (0.500000 , -0.000000) 


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



Init Arrays signal(angle) = cos(angle) + i sin(angle) for freq 6.000000 Hz ... 



Compute DFT  for freq 6.000000 Hz ... 

    Freq[6] = (1.000000 , -0.000000) 


Compute FFT for freq 6.000000 Hz... 

    Freq[6] = (1.000000 , -0.000000) 


    ======================= 



Init Arrays signal(angle) = cos(angle) for freq 7.000000 Hz ... 



Compute DFT  for freq 7.000000 Hz ... 

    Freq[7] = (0.500000 , -0.000000) 


Compute FFT for freq 7.000000 Hz... 

    Freq[7] = (0.500000 , 0.000000) 


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



Init Arrays signal(angle) = cos(angle) + i sin(angle) for freq 7.000000 Hz ... 



Compute DFT  for freq 7.000000 Hz ... 

    Freq[7] = (1.000000 , -0.000000) 


Compute FFT for freq 7.000000 Hz... 

    Freq[7] = (1.000000 , 0.000000) 


    ======================= 


It exist something like a Fast Discret Cosinus Transform, a Fast Discret Sinus Transform and/or a “Fast Discret Cosinus Transform for the real part mixed with a Fast Sinus Transform for the imaginary part” ?

And/or can I use the real part for the left channel and the imaginary part for the right channel on a stereo sound for example ?
(if this is only a story of multiply by two for to have values that I want, it is not really a problem)

[QUOTE=The Little Body;1292046]
The real part is only 0.5f and the imaginary part is 0.0f at the wanted frequency when I use the simple signal(angle) = cos( angle )[/quote]
For a real signal, the DFT has conjugate symmetry: y[N-i]=y[i]*.

Note that eix+e-ix=2·cos(x) and eix-e-ix=2i·sin(x). Also note that eix has period 2π, so e-ix=e2π-ix. So when you take a signal with frequency f and add the image at F-f (where F is the sampling frequency), the imaginary parts cancel and the real part is doubled.

If you’re using this for a spectral display, you just want the magnitude: |z|=sqrt(Re(z)2+Im(z)2). The phase can be ignored.

At first view the thing that I search seem to be the Short Time Fourrier Transform, cf. STFT like explained at https://en.wikipedia.org/wiki/Short-time_Fourier_transform#Discrete-time_STFT

=> this speak of the use of the STFT for to make a spectrogram, so this is relatively certainly what I search :slight_smile: :slight_smile:

Thanks[b] GClements;

If you’re using this for a spectral display, you just want the magnitude: |z|=sqrt(Re(z)2+Im(z)2). The phase can be ignored.

=> I think that is what I search :slight_smile: :slight_smile:
[/b]

I have replaced the print of values computed by the FFT but this is always only the half of 1.0f that I find, cf. 0.5f instead 1.0f


void PrintArrays( int n, bool zeroes )
{
    int i;

    for( i = 0 ; i < n ; i++ )
    {
        if( (fabs(reals[i]) > EPSILON ) || (fabs(images[i]) > EPSILON ) || ( zeroes == true ) )
        {
            printf("    Freq[%d] = (%f , %f) [magnitude=%f] 
", 
                i, reals[i], images[i], sqrt( pow(reals[i], 2.0f) + pow(images[i], 2.0f) ) );
        }
    }
}  

This give this now
(the problem of 0.5f values instead 1.0f is always here :frowning: )


FFT v0.1 by Cyclone 


Init Arrays signal(angle) = cos(angle) for freq 1.000000 Hz ... 



Compute DFT  for freq 1.000000 Hz ... 

    Freq[1] = (0.500000 , -0.000000) [magnitude=0.500000] 


Compute FFT for freq 1.000000 Hz... 

    Freq[1] = (0.500000 , -0.000000) [magnitude=0.500000] 


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



Init Arrays signal(angle) = cos(angle) + i sin(angle) for freq 1.000000 Hz ... 



Compute DFT  for freq 1.000000 Hz ... 

    Freq[1] = (1.000000 , -0.000000) [magnitude=1.000000] 


Compute FFT for freq 1.000000 Hz... 

    Freq[1] = (1.000000 , 0.000000) [magnitude=1.000000] 


    ======================= 



Init Arrays signal(angle) = cos(angle) for freq 2.000000 Hz ... 



Compute DFT  for freq 2.000000 Hz ... 

    Freq[2] = (0.500000 , -0.000000) [magnitude=0.500000] 


Compute FFT for freq 2.000000 Hz... 

    Freq[2] = (0.500000 , 0.000000) [magnitude=0.500000] 


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



Init Arrays signal(angle) = cos(angle) + i sin(angle) for freq 2.000000 Hz ... 



Compute DFT  for freq 2.000000 Hz ... 

    Freq[2] = (1.000000 , -0.000000) [magnitude=1.000000] 


Compute FFT for freq 2.000000 Hz... 

    Freq[2] = (1.000000 , 0.000000) [magnitude=1.000000] 


    ======================= 



Init Arrays signal(angle) = cos(angle) for freq 3.000000 Hz ... 



Compute DFT  for freq 3.000000 Hz ... 

    Freq[3] = (0.500000 , -0.000000) [magnitude=0.500000] 


Compute FFT for freq 3.000000 Hz... 

    Freq[3] = (0.500000 , 0.000000) [magnitude=0.500000] 


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



Init Arrays signal(angle) = cos(angle) + i sin(angle) for freq 3.000000 Hz ... 



Compute DFT  for freq 3.000000 Hz ... 

    Freq[3] = (1.000000 , -0.000000) [magnitude=1.000000] 


Compute FFT for freq 3.000000 Hz... 

    Freq[3] = (1.000000 , -0.000000) [magnitude=1.000000] 


    ======================= 



Init Arrays signal(angle) = cos(angle) for freq 4.000000 Hz ... 



Compute DFT  for freq 4.000000 Hz ... 

    Freq[4] = (0.500000 , -0.000000) [magnitude=0.500000] 


Compute FFT for freq 4.000000 Hz... 

    Freq[4] = (0.500000 , 0.000000) [magnitude=0.500000] 


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



Init Arrays signal(angle) = cos(angle) + i sin(angle) for freq 4.000000 Hz ... 



Compute DFT  for freq 4.000000 Hz ... 

    Freq[4] = (1.000000 , -0.000000) [magnitude=1.000000] 


Compute FFT for freq 4.000000 Hz... 

    Freq[4] = (1.000000 , 0.000000) [magnitude=1.000000] 


    ======================= 



Init Arrays signal(angle) = cos(angle) for freq 5.000000 Hz ... 



Compute DFT  for freq 5.000000 Hz ... 

    Freq[5] = (0.500000 , -0.000000) [magnitude=0.500000] 


Compute FFT for freq 5.000000 Hz... 

    Freq[5] = (0.500000 , 0.000000) [magnitude=0.500000] 


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



Init Arrays signal(angle) = cos(angle) + i sin(angle) for freq 5.000000 Hz ... 



Compute DFT  for freq 5.000000 Hz ... 

    Freq[5] = (1.000000 , -0.000000) [magnitude=1.000000] 


Compute FFT for freq 5.000000 Hz... 

    Freq[5] = (1.000000 , 0.000000) [magnitude=1.000000] 


    ======================= 



Init Arrays signal(angle) = cos(angle) for freq 6.000000 Hz ... 



Compute DFT  for freq 6.000000 Hz ... 

    Freq[6] = (0.500000 , -0.000000) [magnitude=0.500000] 


Compute FFT for freq 6.000000 Hz... 

    Freq[6] = (0.500000 , -0.000000) [magnitude=0.500000] 


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



Init Arrays signal(angle) = cos(angle) + i sin(angle) for freq 6.000000 Hz ... 



Compute DFT  for freq 6.000000 Hz ... 

    Freq[6] = (1.000000 , -0.000000) [magnitude=1.000000] 


Compute FFT for freq 6.000000 Hz... 

    Freq[6] = (1.000000 , -0.000000) [magnitude=1.000000] 


    ======================= 



Init Arrays signal(angle) = cos(angle) for freq 7.000000 Hz ... 



Compute DFT  for freq 7.000000 Hz ... 

    Freq[7] = (0.500000 , -0.000000) [magnitude=0.500000] 


Compute FFT for freq 7.000000 Hz... 

    Freq[7] = (0.500000 , 0.000000) [magnitude=0.500000] 


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



Init Arrays signal(angle) = cos(angle) + i sin(angle) for freq 7.000000 Hz ... 



Compute DFT  for freq 7.000000 Hz ... 

    Freq[7] = (1.000000 , -0.000000) [magnitude=1.000000] 


Compute FFT for freq 7.000000 Hz... 

    Freq[7] = (1.000000 , 0.000000) [magnitude=1.000000] 


    ======================= 



=> it’s normal ?

Note that eix+e-ix=2·cos(x) and eix-e-ix=2i·sin(x). Also note that eix has period 2π, so e-ix=e2π-ix. So when you take a signal with frequency f and add the image at F-f (where F is the sampling frequency), the imaginary parts cancel and the real part is doubled.

=> if I understand good, it’s normal :slight_smile:

Lot of thanks for your help, GClements

This don’t seem to work :frowning:

I have replace the display of FFT coefficients like this :


void PrintArrays( int n, bool zeroes )
{
    int i;

    double realpart;
    double imgpart;

    for( i = 0 ; i < n ; i++ )
    {
        if( (fabs(reals[i]) > EPSILON ) || (fabs(images[i]) > EPSILON ) || ( zeroes == true ) )
        {
            realpart = reals[i] + reals[size-i];
            imgpart  = images[i] + images[size-i];

            printf("    Freq[%d] = (%f , %f) [magnitude=%f] 
", 
                i, realpart, imgpart, sqrt( pow(realpart, 2.0f) + pow(imgpart, 2.0f) ) );
        }
    }
}  

But this give this :


FFT v0.1 by Cyclone 


Init Arrays signal(angle) = cos(angle) for freq 1.000000 Hz ... 



Compute DFT  for freq 1.000000 Hz ... 

    Freq[1] = (1.000000 , -0.000000) [magnitude=1.000000] 


Compute FFT for freq 1.000000 Hz... 

    Freq[1] = (1.000000 , 0.000000) [magnitude=1.000000] 


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



Init Arrays signal(angle) = cos(angle) + i sin(angle) for freq 1.000000 Hz ... 



Compute DFT  for freq 1.000000 Hz ... 

    Freq[1] = (1.000000 , 0.000000) [magnitude=1.000000] 


Compute FFT for freq 1.000000 Hz... 

    Freq[1] = (1.000000 , 0.000000) [magnitude=1.000000] 


    ======================= 



Init Arrays signal(angle) = cos(angle) for freq 2.000000 Hz ... 



Compute DFT  for freq 2.000000 Hz ... 

    Freq[2] = (1.000000 , -0.000000) [magnitude=1.000000] 


Compute FFT for freq 2.000000 Hz... 

    Freq[2] = (1.000000 , 0.000000) [magnitude=1.000000] 


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



Init Arrays signal(angle) = cos(angle) + i sin(angle) for freq 2.000000 Hz ... 



Compute DFT  for freq 2.000000 Hz ... 

    Freq[2] = (1.000000 , -0.000000) [magnitude=1.000000] 


Compute FFT for freq 2.000000 Hz... 

    Freq[2] = (1.000000 , -0.000000) [magnitude=1.000000] 


    ======================= 



Init Arrays signal(angle) = cos(angle) for freq 3.000000 Hz ... 



Compute DFT  for freq 3.000000 Hz ... 

    Freq[3] = (1.000000 , -0.000000) [magnitude=1.000000] 


Compute FFT for freq 3.000000 Hz... 

    Freq[3] = (1.000000 , 0.000000) [magnitude=1.000000] 


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



Init Arrays signal(angle) = cos(angle) + i sin(angle) for freq 3.000000 Hz ... 



Compute DFT  for freq 3.000000 Hz ... 

    Freq[3] = (1.000000 , -0.000000) [magnitude=1.000000] 


Compute FFT for freq 3.000000 Hz... 

    Freq[3] = (1.000000 , -0.000000) [magnitude=1.000000] 


    ======================= 



Init Arrays signal(angle) = cos(angle) for freq 4.000000 Hz ... 



Compute DFT  for freq 4.000000 Hz ... 

    Freq[4] = (1.000000 , -0.000000) [magnitude=1.000000] 


Compute FFT for freq 4.000000 Hz... 

    Freq[4] = (1.000000 , 0.000000) [magnitude=1.000000] 


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



Init Arrays signal(angle) = cos(angle) + i sin(angle) for freq 4.000000 Hz ... 



Compute DFT  for freq 4.000000 Hz ... 

    Freq[4] = (1.000000 , -0.000000) [magnitude=1.000000] 


Compute FFT for freq 4.000000 Hz... 

    Freq[4] = (1.000000 , 0.000000) [magnitude=1.000000] 


    ======================= 



Init Arrays signal(angle) = cos(angle) for freq 5.000000 Hz ... 



Compute DFT  for freq 5.000000 Hz ... 

    Freq[5] = (1.000000 , -0.000000) [magnitude=1.000000] 


Compute FFT for freq 5.000000 Hz... 

    Freq[5] = (1.000000 , 0.000000) [magnitude=1.000000] 


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



Init Arrays signal(angle) = cos(angle) + i sin(angle) for freq 5.000000 Hz ... 



Compute DFT  for freq 5.000000 Hz ... 

    Freq[5] = (1.000000 , -0.000000) [magnitude=1.000000] 


Compute FFT for freq 5.000000 Hz... 

    Freq[5] = (1.000000 , 0.000000) [magnitude=1.000000] 


    ======================= 



Init Arrays signal(angle) = cos(angle) for freq 6.000000 Hz ... 



Compute DFT  for freq 6.000000 Hz ... 

    Freq[6] = (1.000000 , -0.000000) [magnitude=1.000000] 


Compute FFT for freq 6.000000 Hz... 

    Freq[6] = (1.000000 , 0.000000) [magnitude=1.000000] 


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



Init Arrays signal(angle) = cos(angle) + i sin(angle) for freq 6.000000 Hz ... 



Compute DFT  for freq 6.000000 Hz ... 

    Freq[6] = (1.000000 , -0.000000) [magnitude=1.000000] 


Compute FFT for freq 6.000000 Hz... 

    Freq[6] = (1.000000 , -0.000000) [magnitude=1.000000] 


    ======================= 



Init Arrays signal(angle) = cos(angle) for freq 7.000000 Hz ... 



Compute DFT  for freq 7.000000 Hz ... 

    Freq[7] = (1.000000 , -0.000000) [magnitude=1.000000] 


Compute FFT for freq 7.000000 Hz... 

    Freq[7] = (1.000000 , 0.000000) [magnitude=1.000000] 


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



Init Arrays signal(angle) = cos(angle) + i sin(angle) for freq 7.000000 Hz ... 



Compute DFT  for freq 7.000000 Hz ... 

    Freq[7] = (1.000000 , -0.000000) [magnitude=1.000000] 


Compute FFT for freq 7.000000 Hz... 

    Freq[7] = (1.000000 , 0.000000) [magnitude=1.000000] 


    ======================= 

=> I don’t understand why cos(angle) and cos(angle) + i sin(angle) give now exactly the same FFTs :frowning:

==> we “loose” the imaginery part of the original complex signal when we compute the FFT ???

But OK, computed non zeroes coefficients are now good 1.0f instead 0.5f :slight_smile: :slight_smile:

Because adding the coefficients alone simply gives you the real part.

If you add the complete terms, you get:

(a+bi).e-ix
= (a+bi).(cos(-x)+i.sin(-x))
= (a+bi).(cos(x)-i.sin(x))
= (a+bi).cos(x)-(a+bi).i.sin(x)
= a.cos(x)+bi.cos(x)-ai.sin(x)+b.sin(x)
= (a.cos(x)+b.sin(x))-i.(a.sin(x)-b.cos(x))

(a-bi).eix
= (a-bi).(cos(x)+i.sin(x))
= (a-bi).cos(x)+(a-bi).i.sin(x)
= a.cos(x)-bi.cos(x)+ai.sin(x)+b.sin(x)
= (a.cos(x)+b.sin(x))+i.(a.sin(x)-b.cos(x))

Note that these two expressions are conjugates (identical real parts, imaginary parts have the same magnitude but opposite sign).

Adding them together gives:

(a+bi).e-ix + (a-bi).eix
= (a.cos(x)+b.sin(x))-i.(a.sin(x)-b.cos(x)) + (a.cos(x)+b.sin(x))+i.(a.sin(x)-b.cos(x))
= (a.cos(x)+b.sin(x))+(a.cos(x)+b.sin(x)) - i.(a.sin(x)-b.cos(x))+i.(a.sin(x)-b.cos(x))
= 2.(a.cos(x)+b.sin(x)) + 0.i

Note that a.cos(x)+b.sin(x) = r.sin(x+p) where r2=a2+b2 and tan§=a/b.

But if you only want the magnitude, you can just find the magnitude of one of the coefficients and double it.

Essentially: if you graph y=e^ix (one axis for real x, two for complex y) you get a helix. If you graph y=e^-ix you get a helix with the opposite rotation (CW vs CCW). In much the same way that, for real numbers, f(x)=k+x is a translation by k and f(x)=kx is a scaling, for complex numbers f(z)=k+z is a (2D) translation while f(z)=kz is a rotation and scaling. So multiplying a helix by a complex amplitude rotates it, and rotating a helix is the same as translating it (i.e. a phase shift). If you multiply the helix e^-ix by (a+bi) you shift it in one direction; if you multiply the opposing helix e^ix by the conjugate (a-bi) you shift it in the same direction as the first (like turning a right-handed screw clockwise and a left-handed screw counter-clockwise). Regardless of the values a and b, the two always shift by the same amount, and their sum always results in the imaginary parts cancelling out to produce a real sinusoid.

Thanks for yours explanations, I understand now a lot more that before how this work

I have first thinking about using reals parts for to represent the left channel and the imaginary parts for the right channel on a stereo signal but this can’t work as simply as this :frowning:

I have too dream before about the use of something like quaternions for to can handle a quadriphonic channels signal but now I know that is not possible because a more simple stereo signal can not to be handled on a FFT like “simples” complex values :frowning:

But if you only want the magnitude, you can just find the magnitude of one of the coefficients and double it.

It’s now what I think to make because a simple multiply by two is very more speed than a squared root of the sum of two squared value :slight_smile:

Can I delete alls computations on the y array into the FFT( int dir, int m, double *x, double *y ) call for to have a call that is about 2x more speed but with the same x output values because it work only with reals values and imaginary parts are “auto-canceled” because opposites ?

You need to do that anyhow to get the magnitude. If you want the square of the magnitude (which is quite common, as it represents the power), you can omit the square root and just multiply by 4.

[QUOTE=The Little Body;1292061]
Can I delete alls computations on the y array into the FFT( int dir, int m, double *x, double *y ) call for to have a call that is about 2x more speed but with the same x output values because it work only with reals values and imaginary parts are “auto-canceled” because opposites ?[/QUOTE]
No. It is possible to optimise the FFT in the case of a real signal. You still get a complex output, but it only has N/2+1 elements (e.g. 256 real input samples gives you 129 complex amplitude samples).

FWIW, here is some GLSL FFT code I dug up. It isn’t optimised at all.

Compute shader:


#version 430

layout(local_size_x=16) in;
layout(local_size_y=16) in;

layout(rg32f, binding=0) readonly uniform image2D grid_in;
layout(rg32f, binding=1) writeonly uniform image2D grid_out;

uniform int N;
uniform int Ns;
uniform int dir;

const float pi = 3.141592653589793;

vec2 cmul(vec2 a, vec2 b)
{
    return vec2(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x);
}

void main()
{
    ivec2 ij = ivec2(gl_GlobalInvocationID.xy);
    int i = ij.x;
    int j = ij.y;
    int k = Ns / 2;
    int a = (i / Ns) * k;
    int b = i % k;
    int i0 = a + b;
    int i1 = i0 + (N / 2);
    float t = 2*pi*i/Ns;
    vec2 f = vec2(cos(t), sin(t));
    f.y *= dir;
    vec2 x0 = imageLoad(grid_in, ivec2(i0, j)).xy;
    vec2 x1 = imageLoad(grid_in, ivec2(i1, j)).xy;
    imageStore(grid_out, ij, vec4(x0 + cmul(f, x1), 0, 0));
}

Application code (Python):


def do_fft(dir, tex):
    w, h = 512, 512
    gw, gh = 16, 16 # local size

    shader = glCreateShader(GL_COMPUTE_SHADER)
    glShaderSource(shader, fft_x_src)
    glCompileShader(shader)
    if not glGetShaderiv(shader, GL_COMPILE_STATUS):
        raise RuntimeError(glGetShaderInfoLog(shader))

    program = glCreateProgram()
    glAttachShader(program, shader)
    glLinkProgram(program)
    if not glGetProgramiv(program, GL_LINK_STATUS):
        raise RuntimeError(glGetProgramInfoLog(program))

    glUseProgram(program)

    tex1,tex2 = glGenTextures(2)
    glBindTexture(GL_TEXTURE_2D, tex1)
    glTexStorage2D(GL_TEXTURE_2D, 1, GL_RG32F, w, h)
    glBindTexture(GL_TEXTURE_2D, tex2)
    glTexStorage2D(GL_TEXTURE_2D, 1, GL_RG32F, w, h)
    glBindTexture(GL_TEXTURE_2D, 0)

    log_w = int(math.ceil(math.log(w, 2)))

    u_dir = glGetUniformLocation(program, "dir")
    u_N   = glGetUniformLocation(program, "N")
    u_Ns  = glGetUniformLocation(program, "Ns")

    glUniform1i(u_dir, dir)

    for logstep in xrange(log_w):
        src, dst = tex1, tex2
        if logstep == 0:
            src = tex
        step = 2 << logstep
        glBindImageTexture(0, src, 0, False, 0, GL_READ_ONLY,  GL_RG32F)
        glBindImageTexture(1, dst, 0, False, 0, GL_WRITE_ONLY, GL_RG32F)
        glUniform1i(u_N , w)
        glUniform1i(u_Ns, step)
        nx, ny = w / gw, h / gh
        glDispatchCompute(nx, ny, 1)
        tex1, tex2 = tex2, tex1

    glBindImageTexture(0, 0, 0, False, 0, GL_READ_ONLY,  GL_R8)
    glBindImageTexture(1, 0, 0, False, 0, GL_READ_ONLY,  GL_R8)

    glUseProgram(0)

    return dst

Input and output are GL_RG32F textures (R = real, G = imaginary); dir should be -1 for forward FFT, 1 for inverse FFT. Each row of the output is the FFT of the corresponding row of the input. The results agree with NumPy’s FFT to within ~0.05% for most of the range (bear in mind that the GPU is using single precision).