rj-action-library/Runtime/Math/FFT.cs

using Godot;

namespace Rokojori
{
  public class FFT 
  {
    class FloatSineTable 
    {
      public float[] sineValues;

      public FloatSineTable( int numBits ) 
      {
        int len = 1 << numBits;
        sineValues = new float[1 << numBits];

        for ( int i = 0; i < len; i++ ) 
        {
          sineValues[i] = ( float ) Mathf.Sin( ( i * Mathf.Pi * 2.0 ) / len );
        }
      }
    }

    class BitReverseTable 
    {
      public int[] reversedBits;

      public BitReverseTable( int numBits ) 
      {
        reversedBits = new int[1 << numBits];

        for ( int i = 0; i < reversedBits.Length; i++ ) 
        {
          reversedBits[i] = ReverseBits( i, numBits );
        }
      }

      static int ReverseBits( int index, int numBits ) 
      {
        int i = 0;
        int rev = 0;

        for ( i = rev = 0; i < numBits; i++ ) 
        {
          rev = ( rev << 1 ) | ( index & 1 );
          index >>= 1;
        }

        return rev;
      }
    }

    static readonly int MaxSizeLog2 = 16;

    static BitReverseTable[] reverseTables  = new BitReverseTable[ MaxSizeLog2 ];
    static FloatSineTable[] floatSineTables = new FloatSineTable[ MaxSizeLog2 ];
    

    static float[] GetFloatSineTable( int n ) 
    {
      FloatSineTable sineTable = floatSineTables[n];

      if ( sineTable == null ) 
      {
        sineTable = new FloatSineTable( n );

        floatSineTables[n] = sineTable;
      }

      return sineTable.sineValues;
    }
    
    static int[] GetReverseTable( int n ) 
    {
      BitReverseTable reverseTable = reverseTables[n];

      if ( reverseTable == null ) 
      {
        reverseTable = new BitReverseTable( n );
        reverseTables[n] = reverseTable;
      }

      return reverseTable.reversedBits;
    }
    

    static void Transform( int sign, int n, float[] real, float[] imaginary ) 
    {
      float scale = ( sign > 0 ) ? ( 2.0f / n ) : ( 0.5f );

      int numBits = Rokojori.FFT.NumBits( n );

      int[] reverseTable = GetReverseTable( numBits );
      float[] sineTable = GetFloatSineTable( numBits );

      int mask = n - 1;
      int cosineOffset = n / 4; // phase offset between cos and sin

      int i, j;

      for ( i = 0; i < n; i++ ) 
      {
        j = reverseTable[i];

        if ( j >= i ) 
        {
          float tempr = real[j] * scale;
          float tempi = imaginary[j] * scale;
          real[j] = real[i] * scale;
          imaginary[j] = imaginary[i] * scale;
          real[i] = tempr;
          imaginary[i] = tempi;
        }
      }

      int mmax, stride;
      int numerator = sign * n;

      for ( mmax = 1, stride = 2 * mmax; mmax < n; mmax = stride, stride = 2 * mmax ) 
      {
        int phase = 0;
        int phaseIncrement = numerator / ( 2 * mmax );

        for ( int m = 0; m < mmax; ++m ) 
        {
          float wr = sineTable[( phase + cosineOffset ) & mask]; // cosine
          float wi = sineTable[phase];

          for ( i = m; i < n; i += stride ) 
          {
            j = i + mmax;
            float tr = ( wr * real[j] ) - ( wi * imaginary[j] );
            float ti = ( wr * imaginary[j] ) + ( wi * real[j] );
            real[j] = real[i] - tr;
            imaginary[j] = imaginary[i] - ti;
            real[i] += tr;
            imaginary[i] += ti;
          }

          phase = ( phase + phaseIncrement ) & mask;
        }

        mmax = stride;
      }
    }
    
    public static void CalculateMagnitudes( float[] ar, float[] ai, float[] magnitudes ) 
    {
      for ( int i = 0; i < magnitudes.Length; ++i ) 
      {
        magnitudes[i] = ( float ) Mathf.Sqrt( (ar[i] * ar[i] ) + ( ai[i] * ai[i] ));
      }
    }

    static int NumBits( int powerOf2 ) 
    {
      int i;
      
      for ( i = -1; powerOf2 > 0; powerOf2 = powerOf2 >> 1, i++ )
        ;
      return i;
    }

    public static void Forward( int n, float[] real, float[] imaginary ) 
    {
      Transform( 1, n, real, imaginary ); 
    }

    public static void Forward( int n, float[] real ) 
    {
      float[] ai = new float[ real.Length ];

      Forward( n, real, ai );
    }

    public static void Inverse( int n, float[] real, float[] imaginary ) 
    {
      Transform( -1, n, real, imaginary ); 
    }

    public static float FrequencyAtBin( int binIndex, int numBands, float sampleRate )
    {
      return binIndex * sampleRate / ( 2.0f * numBands );
    }

  
  }


}