First you need to figure out which FFT bin your fundamental frequency is in. Say it resides in bin# 10. The harmonics will reside in integer multiples of that bin so the 2nd harmonic will be in bin 20, 3rd in bin 30 and so on. For each of these harmonic bins you need to compute the amplitude. Depending on the window function you used in the FFT you will need to include a small number of bins in the calculation (google spectral leakage if you're interested).
double computeAmpl(double[] spectrum, int windowHalfLen, int peakBin, int harmonic)
{
double sumOfSquares = 0.0;
for (int bin = peakBin-windowHalfLen; bin <= peakBin+windowHalfLen; bin++)
{
sumOfSquares += spectrum[bin] * spectrum[bin];
}
return sqrt(sumOfSquares);
}
As I mentioned the window half length depends on the window. Some common ones are:
- blackman-harris 3 - 3
- blackman-harris 4 - 4
- flat top - 5
- hann - 3