2

I have a Complex[] (from CsCore) which is the result of my FFT.

Complex has a float real and a float imaginary.

From this I calculated the following

  • Frequency: (double)index * sampleRate / FftSize;
  • Amplitude / Magnitude: Math.Sqrt(Math.Pow(real, 2) + Math.Pow(imaginary, 2));
  • Phase: Math.Atan(imaginary / real);

If these are wrong please correct me.

From what I understand, this is the frequency domain information and it allows me to see which frequencies are most common in my sample. Now I want to see the power density over time. Matlab documentation shows examples, But I don't understand it because I don't know Matlab. Could someone explain the Matlab documentation on this subject or help me with a C# implementation?

EDIT:
This answer suggest to simply square the amplitude. Is that correct?

Community
  • 1
  • 1
Gert Kommer
  • 1,163
  • 2
  • 22
  • 49
  • According to [the CsCore `Complex` code](https://github.com/filoe/cscore/blob/master/CSCore/Utils/Complex.cs), it has `float real` and `float imaginary`, not `double`... – JHBonarius Mar 15 '17 at 15:47

2 Answers2

2

Indeed, as I stated in this other answer you could obtain a power spectrum density (PSD) estimate by squaring the amplitudes of the FFT results. This is essentially what the following line from the Matlab documentation you quoted states (up to a scaling factor, which is not significant for most applications requiring only to compare the relative strength of the different frequency components):

psdx = (1/(Fs*N)) * abs(xdft).^2;

As I also mentioned in my other answer, and is also described in Matlab documentation, you could get a better PSD estimate by multiplying your signal by a window function before taking the FFT, and averaging the squared-magnitude of multiple FFT results.

Note: for the phase you would be better served with Math.Atan2(imaginary, real) (see Math.Atan2 on MSDN) which covers the enter [-pi,pi] range (instead of Math.Atan() which only covers [-pi/2,pi/2]).

Community
  • 1
  • 1
SleuthEye
  • 14,379
  • 2
  • 32
  • 61
  • I'm not really sure how to use the windowing. I understand the formula of the hann(`0.50f * (1 - (float)Math.Cos((2 * Math.PI * n) / N - 1))`) but I have to use the window before the fft on my wave? – Gert Kommer Mar 15 '17 at 15:16
  • 1
    If you just use the [`Complex.Phase` property](https://msdn.microsoft.com/en-us/library/system.numerics.complex.phase(v=vs.110).aspx) it already uses `Math.Atan2(b, a)` for a complex number `a + bi`. – JHBonarius Mar 15 '17 at 15:36
  • @GertKommer yes, if you have a chunk of N samples from your wave, you'd multiply element wise by the window function, then take the FFT of the result. – SleuthEye Mar 15 '17 at 15:42
  • I have a `float[] buffer1` which contains 10ms of data. I just have to create a new `float[] buffer2` that holds the values of `buffer1 * window`. So if `buffer1[0] == 5` and output of the window function is `2` than `buffer2[0] = 10` and so on – Gert Kommer Mar 15 '17 at 16:06
1

Firstly Math.Sqrt(Math.Pow(real, 2) + Math.Pow(imaginary, 2)); is already implemented as the Complex.Magnitude property. Or you can use the Complex.Abs method.

In addition to what SleuthEye said, I did some measurements on function implementation.

Because I did not trust the Math.Pow(x,2) function I implemented:

private static double Square(double value)
{
    return value * value;
}

However, it turns out that C# already optimized Math.Pow(x,2), so it's fast enough. But anyhow: next I compared three implementations

  1. Square(testData[idx].Real) + Square(testData[idx].Imaginary);
  2. Square(testData[idx].Magnitude);
  3. Square(Complex.Abs(testData[idx]));

My (average) results were (for 10,000,000 complex elements):

  1. 45 ms
  2. 220 ms
  3. 211 ms

So it seems the Magnitude property and Abs method use a square root internally, which takes a lot of cycles to process. But for the PSD you don't need that.

JHBonarius
  • 10,824
  • 3
  • 22
  • 41
  • Im not using the standard Complex class, but the Complex class of CsCore. Sadly that one doesn't have the magnitude and phase. As for the rest of your answer, I'm not really sure what to do with it. – Gert Kommer Mar 15 '17 at 15:37
  • Sometimes its worth it to convert one type to another to be able to use more functions ;) But what I am trying to say: first determining the Absolute/Magnitude (I see: in CsCore it is called `Complex.Value`) -which has a square root- and then squaring the result is very inefficient. – JHBonarius Mar 15 '17 at 15:44
  • Well I'm ashamed I didn't see the `Complex.Value`. For the PSD i do have to square over the magnitude and then the `Complex.Value` or my own writen function is faster than the Sys numerics? – Gert Kommer Mar 15 '17 at 16:03
  • `Math.Pow(a[i].Real, 2) + Math.Pow(a[i].Imaginary, 2)` is faster then `Math.Pow(a[i].Value, 2)`. It is a good exercise to write a comparison test yourself :) – JHBonarius Mar 15 '17 at 16:14
  • 1
    Well the `a[i].Value` also has a sqrt and the other one does not, so that it is faster is not really a suprise. Good suggestion about comparing different methods myself. Honestly I didn't do it because I was more interested in how the implementation worked and in a later phase I'm going to worry about performance and start optimizing. But still thx for the tip:) – Gert Kommer Mar 15 '17 at 16:24