using discrete fourier transform to observe pure sinusoidal signal

Mar 6, 2013 at 11:33 PM
hello,

I want to use discrete fourier tranform with C#. so mathnet seemed like a perfect candidate for the job.

During a 5 seconds capture, i can capture about 40~57 values (sometimes more, sometimes less).

I know the captured data is expected to be a sinusoid. The period is 1 sec so there should be about 5 periods in the captured data. I want to give a quality score to the measure (how far from the model) :
0 = bad.
1 = perfect score

my problem is that the acquisition is slow/unreliable and during the same time frame (5sec) i get more or less samples between runs.

I believe that in order to give a "quality score" to such a feed i should focus on :
  • make sure the maximum values are located in the correct place in Complex[] (= observed frequency is ok)
  • make sure the ratio between the max and the biggest other value decent (= signal to noise ratio is ok)
However, the use of Complex[] samples in mathnet confuses me + FourierOptions confuses me too + i am unsure i should use DTF in this you-never-have-same-amount-of-samples situation + i have not made any math in over 10 years. i need help :)

My questions
  • is the approach correct ? does it make sense?
  • is using Complex(mydata, 0) ok for each input point ? (ignoring imaginary value)
  • am i right to ignore phase and use only magnitude in the Complex[] result ?
thanks

alex

class Program
{
    static void Main(string[] args)
    {
        double[] scores = new double[10];
        scores[0] = fft_test(54, 5.00, null);
        scores[1] = fft_test(47, 5.00, null);
        scores[2] = fft_test(41, 5.00, new[] { 31 });
        scores[3] = fft_test(43, 5.00, new[] { 25 });
        scores[4] = fft_test(55, 5.00, new[] { 12, 22 });
        scores[5] = fft_test(43, 5.00, new[] { 12, 22 });
        scores[6] = fft_test(54, 5.00, new[] { 12, 35, 47 });
        scores[7] = fft_test(47, 5.00, new[] { 12, 35, 47 });

    }

    static double fft_test(int numPoints, double numPeriods, int[] stupidIndexes)
    {
        Random rnd = new Random();
        double[] x_values = new double[numPoints];
        double[] y_values = new double[numPoints];
        //double[] magnitudes = new double[numPoints];
        double start = 0;
        double end = numPeriods * 2 * Math.PI;
        Complex[] dftme = new Complex[numPoints];
        for (int i = 0; i < numPoints ; i++)
        {
            x_values[i] = (double)(i * (end - start)) / numPoints; 
            y_values[i] = (stupidIndexes != null && stupidIndexes.Contains(i) ?  Math.Sin(rnd.NextDouble() *2 * Math.PI)  :   Math.Sin(x_values[i]));
            dftme[i] = new Complex(y_values[i], 0);
        }

        //MathNet.Numerics.IntegralTransforms.Transform.FourierForward(fftme, MathNet.Numerics.IntegralTransforms.FourierOptions.NoScaling); // is this the right fonction ? or do i have to instantiace DiscreteFourierTransform ?
        MathNet.Numerics.IntegralTransforms.Algorithms.DiscreteFourierTransform dft = new MathNet.Numerics.IntegralTransforms.Algorithms.DiscreteFourierTransform();
        dft.BluesteinForward(dftme, MathNet.Numerics.IntegralTransforms.FourierOptions.NoScaling); // what is this FourierOption thing ? i just want the "usual" one

        string res = "";
        for (int i = 0; i < numPoints; i++) { res += y_values[i] + ";" + dftme[i].Real + ";" + dftme[i].Imaginary + ";" + dftme[i].Magnitude + ";" + dftme[i].Phase + "|"; }
        //for (int i = 0; i < numPoints; i++) { magnitudes[i] = dftme[i].Magnitude; }
        //Array.Sort(magnitudes);
        //Array.Reverse(magnitudes);
        int[] topIndices = getPeakIndices(dftme, 5);

        int[] expectedIndices = { (int)numPeriods, (int)(numPoints - numPeriods) };

        double freqScore = (expectedIndices.Contains(topIndices[0]) && expectedIndices.Contains(topIndices[1]) ? 1.0 : 0.0);
        double signalToNoiseScore = 1 - (dftme[topIndices[3]].Magnitude / dftme[topIndices[0]].Magnitude);


        return freqScore * signalToNoiseScore;

    }

    static int[] getPeakIndices(Complex[] dtfme, int numPoints)
    {
        int[] results = new int[numPoints];
        double currentMax = double.MinValue;
        double previousMax = double.MaxValue;
        int currentIndex = -1;
        // this o(n²) smells like bug. it is too late at night for clear thinking
        for (int k = 0; k < numPoints; k++)
        {
            for (int i = 0; i < dtfme.Length; i++)
            {
                if ( currentMax < dtfme[i].Magnitude && dtfme[i].Magnitude < previousMax)
                {
                    currentMax = dtfme[i].Magnitude;
                    currentIndex = i;
                }
            }
            results[k] = currentIndex;
            previousMax = currentMax;
            currentMax = double.MinValue;
        }
        return results;
    }



}




Mar 11, 2013 at 12:29 PM
Edited Mar 11, 2013 at 12:54 PM
  • is the approach correct ? does it make sense?
The approach makes sense, but I think the implementation is flawed (your expectedIndices definition and the following code). Make sure you understand where the frequency peaks should be located in your FFT output (see your next question). Also consider spectral leakage.
  • is using Complex(mydata, 0) ok for each input point ? (ignoring imaginary value)
Yes, but since you're doing a complex DFT on real data, the output frequencies are symmetrical around n/2, so after your FFT you should only consider the first n/2 values for further processing.
  • am i right to ignore phase and use only magnitude in the Complex[] result ?
Yes.


Additionally, your getPeakIndices method isn't correct. Why not just use Array.Sort:
class MagnitudeComparer : IComparer<Complex>
{
    public int Compare(Complex a, Complex b)
    {
    return b.Magnitude.CompareTo(a.Magnitude);
    }
}

static int[] getPeakIndices(Complex[] dtfme, int numPoints)
{
    var indices = Enumerable.Range(0, dtfme.Length).ToArray();
    
    // Considering only the first half of the FFT data.
    Array.Sort(dtfme, indices, 0, dtfme.Length / 2, new MagnitudeComparer());

    return indices.Take(numPoints).ToArray();
}
Mar 11, 2013 at 4:37 PM
Here's a version of getPeakIndices which will work. the above version will also sort the complex values, so the indices won't point to the right location anymore.
static int[] getPeakIndices(Complex[] dtfme, int numPoints)
{
    int n = dtfme.Length / 2;

    var magnitudes = dtfme.Take(n).Select((z) => z.Magnitude).ToArray();
    var indices = Enumerable.Range(0, n).ToArray();

    Array.Sort(magnitudes, indices);

    return indices.Reverse().Take(numPoints).ToArray();
}