1

I am attempting to use Lomont FFT in order to return complex numbers to build a spectrogram / spectral density chart using c#.

I am having trouble understanding how to return values from the class. Here is the code I have put together thus far which appears to be working.

    int read = 0;
    Double[] data;
    byte[] buffer = new byte[1024];

    FileStream wave = new FileStream(args[0], FileMode.Open, FileAccess.Read);  
    read = wave.Read(buffer, 0, 44);                                            
    read = wave.Read(buffer, 0, 1024);                                          
    data = new Double[read];                                                    


    for (int i = 0; i < read; i+=2)
    {
        data[i] = BitConverter.ToInt16(buffer, i) / 32768.0;
        Console.WriteLine(data[i]);
    }
    LomontFFT LFFT = new LomontFFT();
    LFFT.FFT(data, true);

What I am not clear on is, how to return/access the values from Lomont FFT implementation back into my application (console)?

Being pretty new to c# development, I'm thinking I am perhaps missing a fundamental aspect of understanding regarding how to retrieve processed values from the instance of the Lomont Class, or perhaps even calling it incorrectly.

Console.WriteLine(LFFT.A); // Returns 0
Console.WriteLine(LFFT.B); // Returns 1

I have been searching for a code snippet or explanation of how to do this, but so far have come up with nothing that I understand or explains this particular aspect of the issue I am facing. Any guidance would be greatly appreciated.

A subset of the results held in data array noted in the code above can be found below and based on my current understanding, appear to be valid:

0.00531005859375
0.0238037109375
0.041473388671875
0.0576171875
0.07183837890625
0.083465576171875
0.092193603515625
0.097625732421875
0.099639892578125
0.098114013671875
0.0931396484375
0.0848388671875
0.07354736328125
0.05963134765625
0.043609619140625
0.026031494140625
0.007476806640625
-0.011260986328125
-0.0296630859375
-0.047027587890625
-0.062713623046875
-0.076141357421875
-0.086883544921875
-0.09454345703125
-0.098785400390625
-0.0994873046875
-0.0966796875
-0.090362548828125
-0.080810546875
-0.06842041015625
-0.05352783203125
-0.036712646484375
-0.0185546875

What am I actually attempting to do? (perspective)

I am looking to load a wave file into a console application and return a spectrogram/spectral density chart/image as a jpg/png for further processing.

The wave files I am reading are mono in format


UPDATE 1

I Receive slightly different results depending on which FFT is used.

Using RealFFT

    for (int i = 0; i < read; i+=2)
    {
        data[i] = BitConverter.ToInt16(buffer, i) / 32768.0;
        //Console.WriteLine(data[i]);
    }

    LomontFFT LFFT = new LomontFFT();
    LFFT.RealFFT(data, true);

    for (int i = 0; i < buffer.Length / 2; i++)
    {
        System.Console.WriteLine("{0}",
          Math.Sqrt(data[2 * i] * data[2 * i] + data[2 * i + 1] * data[2 * i + 1]));
    }

Partial Result of RealFFT

0.314566983321381
0.625242818210924
0.30314888696868
0.118468857708093
0.0587697011760449
0.0369034115568654
0.0265842582236275
0.0207195964060356
0.0169601273233317
0.0143745438577886
0.012528799609089
0.0111831275153128
0.0102313284519146
0.00960198279358434
0.00920236001619566

Using FFT

    for (int i = 0; i < read; i+=2)
    {
        data[i] = BitConverter.ToInt16(buffer, i) / 32768.0;
        //Console.WriteLine(data[i]);
    }

    double[] bufferB = new double[2 * data.Length];

    for (int i = 0; i < data.Length; i++)
    {
        bufferB[2 * i] = data[i];
        bufferB[2 * i + 1] = 0;
    }

    LomontFFT LFFT = new LomontFFT();
    LFFT.FFT(bufferB, true);

    for (int i = 0; i < bufferB.Length / 2; i++)
    {
        System.Console.WriteLine("{0}",
          Math.Sqrt(bufferB[2 * i] * bufferB[2 * i] + bufferB[2 * i + 1] * bufferB[2 * i + 1]));
    }

Partial Result of FFT:

0.31456698332138
0.625242818210923
0.303148886968679
0.118468857708092
0.0587697011760447
0.0369034115568653
0.0265842582236274
0.0207195964060355
0.0169601273233317
0.0143745438577886
0.012528799609089
0.0111831275153127
0.0102313284519146
0.00960198279358439
0.00920236001619564
Majickal
  • 176
  • 2
  • 16

1 Answers1

0

Looking at the LomontFFT.FFT documentation:

Compute the forward or inverse Fourier Transform of data, with data containing complex valued data as alternating real and imaginary parts. The length must be a power of 2. The data is modified in place.

This tells us a few things. First the function is expecting complex-valued data whereas your data is real. A quick fix for this is to create another buffer of twice the size and setting all the imaginary parts to 0:

double[] buffer = new double[2*data.Length];
for (int i=0; i<data.Length; i++)
{
  buffer[2*i] = data[i];
  buffer[2*i+1] = 0;
}

The documentation also tells us that the computation is done in place. That means that after the call to FFT returns, the input array is replaced with the computed result. You could thus print the spectrum with:

LomontFFT LFFT = new LomontFFT();
LFFT.FFT(buffer, true);

for (int i = 0; i < buffer.Length/2; i++)
{
  System.Console.WriteLine("{0}",
    Math.Sqrt(buffer[2*i]*buffer[2*i]+buffer[2*i+1]*buffer[2*i+1]));
}

Note since your input data is real valued you could also use LomontFFT.RealFFT. In that case, given a slightly different packing rule, you would obtain the FFT results using:

LomontFFT LFFT = new LomontFFT();
LFFT.RealFFT(data, true);

System.Console.WriteLine("{0}", Math.Abs(data[0]);
for (int i = 1; i < data.Length/2; i++)
{
  System.Console.WriteLine("{0}",
    Math.Sqrt(data[2*i]*data[2*i]+data[2*i+1]*data[2*i+1]));
}
System.Console.WriteLine("{0}", Math.Abs(data[1]);

This would give you the non-redundant lower half of the spectrum (Unlike LomontFFT.FFT which provides the entire spectrum). Also, numerical differences on the order of double precision (around 1e-16 times the spectrum peak value) with respect to LomontFFT.FFT can be expected.

SleuthEye
  • 14,379
  • 2
  • 32
  • 61
  • Thanks for the explanation, what you have pointed out makes good sense to me. I appear to be receiving slightly different results depending on which FFT algorithm I use, RealFFT or FFT. I am not sure how significant this actually is though? Updated 1 in the OP describes this. @SleuthEye – Majickal Sep 07 '16 at 03:58
  • @SleutheEye, so to clarify here, if I can only expect nyquist and below results from RealFFT, why is it we are dividing the `data.Length` by 2 `data.Length/2` in the for loop? Is this due to imaginary numbers being returned from RealFFT and the halving of `data.Length` used as a means to filter all items above the nyquist? – Majickal Sep 11 '16 at 11:43
  • 1
    `data.Length/2` for the loop is due to complex numbers being made of an imaginary and real parts (ie. 2 `double` values per complex number). The fact that `RealFFT` only return the lower half of the spectrum implies that `data.Length == buffer.Length/2` (for a corresponding equivalent `buffer` used with `FFT` call). So, if you have 1024 real values, you'd have `data.Length==1024`, `buffer.Length==2048` and you'd execute the `RealFFT` snippet's loop body 511 times (we skip index 0, which contains 2 purely real component at 0Hz and Nyquist frequencies respectively). – SleuthEye Sep 11 '16 at 17:09