
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 youneverhavesameamountofsamples 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();
}



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();
}

