Calculate THD using FHT library

Hello, I’m trying to create a simple audio spectrum analyzer. I’m using LM386 for amplify the signal from the microphone. I’ve done this using the following code and then using processing to visualize the spectrum.

   #define LOG_OUT 1 // use the log output function
    
   #define FHT_N 256 // set to 256 point fht
   
   #include <FHT.h> // include the library

   
   void setup() {
     Serial.begin(115200); // use the serial port
     TIMSK0 = 0; // turn off timer0 for lower jitter
     ADCSRA = 0xe5; // set the adc to free running mode
     ADMUX = 0x40; // use adc0
     DIDR0 = 0x01; // turn off the digital input for adc0
   }
   
   void loop() {
     while(1) { // reduces jitter
       cli();  // UDRE interrupt slows this way down on arduino1.0
       for (int i = 0 ; i < FHT_N ; i++) { // save 256 samples
         while(!(ADCSRA & 0x10)); // wait for adc to be ready
         ADCSRA = 0xf5; // restart adc
         byte m = ADCL; // fetch adc data
         byte j = ADCH;
         int k = (j << 8) | m; // form into an int
         k -= 0x0200; // form into a signed int
         k <<= 6; // form into a 16b signed int
         fht_input[i] = k; // put real data into bins
       }
     


       fht_window(); // window the data for better frequency response
       fht_reorder(); // reorder the data before doing the fht
       fht_run(); // process the data in the fht
       fht_mag_log(); // take the output of the fht
       //fht_mag_octave();
       sei();
       //Serial.write(255); // send a start byte
       //Serial.write(fht_log_out, FHT_N/2); // send out the data
     
   
       for(int i=0;i<FHT_N/2;i++){
        Serial.print("Harmonic");
        Serial.print(i);
        Serial.print(" :");
        Serial.println(fht_log_out[i]);

       }
     
     }
   }

My goal is to measure the THD (total harmonic distortion) but I don’t understand the unit of measure of the harmonics output from the code. I need the RMS value of the harmonics and then calculate the THD according the equation in the following image.

Thanks

I've never used FFT or FHT, but I'll be surprised if you get anything useful... I think you'll get "leakage" between the bins and I don't think you've got enough resolution to see the real distortion, unless the distortion is really bad.

What do the raw results look like? With a "pure" sine wave are you seeing a strong fundamental with nothing (or almost nothing) in the other bins?

I assume the FHT results are approximately an average so that should be "close enough" to RMS or you can multiply the average by 1.11 to get RMS.*

Do you know the frequency, is it fixed, and is it a sine wave? If so, there should be a bin containing the fundamental and from there you can determine which bins contain the harmonics. Square the bin values, sum the results, and take the square root. Then divide by the value from the fundamental bin.

Or, you may want to calculate the RMS for all of the bins which would be the total distortion, not just the harmonic distortion.

  • For a sine wave RMS = 0.707 x peak and average = 0.637 x peak.

What do the raw results look like? With a "pure" sine wave are you seeing a strong fundamental with nothing (or almost nothing) in the other bins?

The raw data are :

Harmonic0 :189
Harmonic1 :173
Harmonic2 :74
Harmonic3 :30
Harmonic4 :46
Harmonic5 :33
Harmonic6 :27
Harmonic7 :56
Harmonic8 :19
Harmonic9 :30
Harmonic10 :39
Harmonic11 :19
Harmonic12 :19
Harmonic13 :33
Harmonic14 :33
Harmonic15 :19
Harmonic16 :8
Harmonic17 :46
Harmonic18 :24
Harmonic19 :16
Harmonic20 :30
Harmonic21 :19
Harmonic22 :30
Harmonic23 :0
Harmonic24 :30
Harmonic25 :19
Harmonic26 :8
Harmonic27 :19
Harmonic28 :8
Harmonic29 :30
Harmonic30 :0
Harmonic31 :37
Harmonic32 :0
Harmonic33 :25
Harmonic34 :19
Harmonic35 :0
Harmonic36 :19
Harmonic37 :8
Harmonic38 :30
Harmonic39 :0
Harmonic40 :27
Harmonic41 :8
Harmonic42 :8
Harmonic43 :8
Harmonic44 :27
Harmonic45 :0
Harmonic46 :30
Harmonic47 :0
Harmonic48 :19
Harmonic49 :8
Harmonic50 :16
Harmonic51 :27
Harmonic52 :27
Harmonic53 :0
Harmonic54 :27
Harmonic55 :19
Harmonic56 :0
Harmonic57 :0
Harmonic58 :16
Harmonic59 :24
Harmonic60 :0
Harmonic61 :19
Harmonic62 :0
Harmonic63 :8
Harmonic64 :8
Harmonic65 :27
Harmonic66 :0
Harmonic67 :19
Harmonic68 :0
Harmonic69 :24
Harmonic70 :16
Harmonic71 :30
Harmonic72 :19
Harmonic73 :0
Harmonic74 :27
Harmonic75 :0
Harmonic76 :8
Harmonic77 :8
Harmonic78 :0
Harmonic79 :27
Harmonic80 :0
Harmonic81 :19
Harmonic82 :33
Harmonic83 :0
Harmonic84 :27
Harmonic85 :8
Harmonic86 :8
Harmonic87 :8
Harmonic88 :8
Harmonic89 :16
Harmonic90 :19
Harmonic91 :8
Harmonic92 :0
Harmonic93 :19
Harmonic94 :0
Harmonic95 :19
Harmonic96 :0
Harmonic97 :0
Harmonic98 :16
Harmonic99 :0
Harmonic100 :8
Harmonic101 :0
Harmonic102 :27
Harmonic103 :0
Harmonic104 :19
Harmonic105 :0
Harmonic106 :19
Harmonic107 :0
Harmonic108 :19
Harmonic109 :16
Harmonic110 :0
Harmonic111 :16
Harmonic112 :8
Harmonic113 :0
Harmonic114 :30
Harmonic115 :8
Harmonic116 :19
Harmonic117 :0
Harmonic118 :19
Harmonic119 :19
Harmonic120 :0
Harmonic121 :0
Harmonic122 :8
Harmonic123 :8
Harmonic124 :16
Harmonic125 :19
Harmonic126 :0
Harmonic127 :8

This because the generic element of the output array fht_log_out is uint8_t. So I have 255 possible levels (8 bit).

From the library (FHTFunctions - Open Music Labs Wiki) I found this :

fht_mag_log() - This gives the magnitude of each bin in the FHT. It sums the squares of the imaginary and real, and then takes the square root, and then takes the log base 2 of that value. Therefore, the output is compressed in a logarithmic fashion, and is essentially in decibels (times a scaling factor). It takes no variables, and returns no variables. It uses a lookup table to calculate the log of the square root, and scales the output over the full 8b range {the equation is 16*(log2((img2 + real2)1/2))}. It is only an 8b value, and the values are taken from fht_input, and returned in fht_log_out. The output values are in sequential order of FHT frequency bins, and there are only N/2 total bins, as the second half of the FFT result is redundant.

So my question is how can I obtain the RMS values of the harmonic from the raw output data show before?

Do you know the frequency, is it fixed, and is it a sine wave?

The frequency is variable, it depends from voice through the microphone.

Or, you may want to calculate the RMS for all of the bins which would be the total distortion, not just the harmonic distortion.

I want calculate THD (total harmonic distortion ) . So I need the RMS value of the all harmonics (128 samples in my case).

The frequency is variable, it depends from voice through the microphone.

Look at the equation in the image you linked to in your first post. The frequency of your voice varies continuously, so how will you measure Vfund_rms and how will you distinguish between a frequency that is part of your voice and one that is a harmonic of your voice?

Harmonic0 :189
Harmonic1 :173
Harmonic2 :74
Harmonic3 :30

Those aren't harmonics. They are FHT bin numbers whose width, in Hz, is determined by the sampling frequency and the number of sampling points in the FHT input. So, first you need to know the sampling frequency with some reasonable degree of accuracy, otherwise you can't make any meaningful measurement of THD. When I played around with the FHT, I got a sampling rate from that code of around 9kHz. This means that the highest frequency in the FHT output is around 4.5kHz.

to measure the THD (total harmonic distortion)

The THD of what? The microphone? The LM386?

Pete

Look at the equation in the image you linked to in your first post. The frequency of your voice varies continuously, so how will you measure Vfund_rms and how will you distinguish between a frequency that is part of your voice and one that is a harmonic of your voice?

I think to calculate the THD in real time,instantly. I know the Vfund_rms is the first harmonic in the output array, and the other elements are the other harmonics to put in the equation.

So, first you need to know the sampling frequency with some reasonable degree of accuracy, otherwise you can't make any meaningful measurement of THD. When I played around with the FHT, I got a sampling rate from that code of around 9kHz. This means that the highest frequency in the FHT output is around 4.5kHz.

My sampling rate is 38.4 Khz so the highest frequency in the FHT output is around 19,2 Khz.

The THD of what? The microphone? The LM386?

I want to calculate the THD of the LM386, but I think it is more difficult. So I decided to calculate the THD of the circuit (LM386+microphone). But the problem is the same how can I obtain RMS value harmonics from FHT bin numbers?

Thanks

I know the Vfund_rms is the first harmonic in the output array

That's all very well, but what is the frequency of the first harmonic of your voice and where is it in the FHT output?

My sampling rate is 38.4 Khz

How do you know that? You would be best to measure the frequency to be sure.

But the problem is the same how can I obtain RMS value harmonics from FHT bin numbers?

The question is the same. What is the frequency of the first harmonic of your voice?
Actually, your (or anyone else's) voice is useless as a test signal for THD. You need to generate a pure sine wave signal of a known frequency into the microphone. Then you can ask about bin numbers.

If the number of samples in the FHT is N and the sampling frequency is f (Hz) then the width of each bin is f/N Hz. The bin numbers are indexed starting at zero where the zero'th bin is the DC component of the signal. Bin number n corresponds to frequency n*f/N Hz. You know the size of the FHT (in your code it is 256). If your sampling rate really is 38.4kHz then each bin in the FHT output is 150Hz wide. This brings up another complication. Is the resolution of the bin frequencies sufficient to determine a meaningful THD?
If you know the frequency of the sine wave test signal, you can easily calculate its bin number and the bin numbers of the harmonics.

Pete

How do you know that? You would be best to measure the frequency to be sure.

I follow this guide Arduino sound level meter and spectrum analyzer | Arik Yavilevich's blog and I have set the ADC in free running and I set the ADC speed using the ADCSRA register ( Each sample takes the ADC about 13 clock cycles to get processed. So by default we get 16Mhz/128/13=9846Hz sampling. If we want to sample at double the rate we can change the divider to be 64 instead.Set divider to 32 which equals a sampling rate of 16Mhz/32/13~=38Khz.

To be sure I have done a test using http://www.szynalski.com/tone-generator/ . I generate a sinusoidal signal with 17 Khz(see the image ). If I want to measure the sampling frequency, how can I do ?

If you know the frequency of the sine wave test signal, you can easily calculate its bin number and the bin numbers of the harmonics.

Ok. Suppose I generate a sine wave with the generator tone. The signal go into the microphone , but now I can use FHT bin to measure THD?

Thanks a lot for help (sorry If I say something wrong)

17kHz isn't a useful frequency for a test of FHT. It's second harmonic (and higher harmonics) is outside the range of the FHT. If the sampling frequency is 38kHz, the highest frequency you can see in the FHT is half that, 19kHz.
You should use a test tone of about 1kHz, for example.

If the test tone really was exactly 17kHz, the output in your FHT_reply.jpg image suggests that the sampling frequency is lower than 38kHz. If the sampling rate is 16MHz/32/13 and the number of samples is 256, the bin width is 150.24Hz. A frequency of 17kHz should fall in bin 17000/150.24 = 113. But your image seems to show peaks in bins 118/119 which means the sampling rate is more like 36.8Hz.
The image also shows up another problem. If the 17kHz tone is smeared across at least 4 bins, as appears to be the case, you aren't going to be able to reliably measure the amplitude of the fundamental frequency, and the harmonics will be even more difficult.

I measured the sampling rate by adding a line of code which toggled the builtin LED on and off and then connected my multimeter to pin 13 - my multimeter has a frequency counter setting.

Pete

My sampling rate is 38.4 Khz so the highest frequency in the FHT output is around 19,2 Khz.

Section 28.1 of the [u]ATmega Datasheet[/u] says you can only get accurate (10-bit) resolution up to a 15kHz sample rate. (7.5kHz audio).

I want to calculate the THD of the LM386, but I think it is more difficult.

In that case use a direct connection (no microphone). The speaker & microphone are going to generate far more distortion than the amplifier, unless there is something terribly wrong with the amplifier.

Your raw data seems to have LOTs of "distortion". But, I doubt that's real... I'm pretty sure it's FHT artifacts. Except for the bin containing the fundamental, all of the other readings should be nearly zero (or zero).

The LM386 datasheet says 0.2% typical THD, which is too low to measure without sophisticated equipment. Although the distortion could be worse if the circuit is built wrong, the chip is defective, or if you are driving it into clipping.

Distortion analyzers are not easy to build and they are expensive. Normally, you'd use a very-very good notch filter to remove the fundamental, then measure what remains. [u]Wikipedia- THD analyzer[/u]

I'd suggest calculating the power spectrum, and then summing the distortion harmonics linearly from that.
From an FFT output you just need to calculate the magnitude-squared for each bin for its power.

The window function used in the FFT will attentuate the signal by a constant factor that you'd need to
correct for to get accurate absolute values, but for THD we are only interested in ratio of power from
fundamental to the rest, so that can be ignored.

You'd have to lookup how to calculate power spectrum from the Hadamard transform, if its even possible...

Another way to calculate THD doesn't require any transform, you estimate the frequency, then correlate
the signal samples with a model sine wave, refining the frequency, phase and amplitude estimates till
it converges on the best match, then you subtract the model from the sampled wave and compare rms
values for signal and residuals directly - this gives THD+noise in one go. Its O(N) rather than O(N logN).

Its a bit more code as you have to do an initial frequency estimate (zero crossings), then have some
correlation code, a sine-wave generator (cordic?), arctangent to determine phase, and code to calculate rms,
so the FFT might be quicker to get working.

Other methods include directly mimicing analog distortion meters by using a tunable digital notch IIR filter
and adjusting it for a null. (this is done in real time and reduces the storage needed to just the few samples needed for the filter).