[ASK] Compute THD and SNR thru teensy 3.2 board only

Status
Not open for further replies.
A

alansoong97259

Member
Hi all,

Currently, I'm using Teensy 3.2 to develop some simple tester to measure freq and Amp(RMS) for speaker and microphone DUT.
And I already completed and verified the code fro the freq and amp.
Now, I wish to have additional functional test with teensy 3.2 on THD and SNR computation.
Do any audio experts out there have this knowledge or experience can share with me?

Thanks in advance.

br
Alan Soong
 
If you know the frequency, maybe a FFT could get you close? If the frequency isn't exactly on a FFT bin, of course you'll need to use one of the windows which causes some smearing. So perhaps take the 2 or 3 closest bins and assume that's all due to the frequency you're measuring? Then multiple the bin number by integers and get the numbers from those groups of bins to find out how much energy is in all the harmonics. If you know the signal was supposed to be just 1 pure sine wave, and if you're using a window that prevents spectral leakage, maybe just add up all the bins not close to the frequency and call that the total noise? Not perfect, but maybe good enough?
 
Hi Paul, Thank you for your support and suggestion.
For my freq measurement, I'm using pulse count with interrupt method to measure.
Please see below my code for the freq and amp detection.

Code: 
// period of pulse accumulation and serial output, milliseconds
#define MainPeriod 128
#define f_input_pin A2   // input pin for pulseIn aka pin 16

long previousMillis = 0; // will store last time of the cycle end
volatile unsigned long duration=0; // accumulates pulse width
volatile unsigned int pulsecount=0;
volatile unsigned long previousMicros=0;

const int dc_offset = 512;          // DC Offset to be removed from incoming signal 
const int numberOfSamples = 128;    // Number of samples to read at a time
const float aref_voltage = 3.3;     // Reference voltage of the Arduino ADC
const float alpha = 0.95;           // smoothing value

int sample;           // individual reading
long signal;     
long averageReading;
long RMS;             // root mean squared calculation

float db;
float dbspl;
long smoothedValue = 0;

void setup()
{
  pinMode(f_input_pin,INPUT);
  Serial.begin(19200);
  attachInterrupt(16, myinthandler, RISING);
  
}

void loop()
{
  if(Serial.available())
  {
      String ch;
      ch = Serial.readString();
      ch.trim();
      if(ch=="start"||ch=="START")
      {
        tone_test();
      }
  }
} 

void tone_test()
{
  duration = 0; // clear counters
  pulsecount = 0;
  tone(8,1000,1000); // Create 1kHz tone for 1 sec on pin 8
  previousMicros = micros();
  unsigned long currentMillis = millis();
  previousMillis = currentMillis;
  unsigned long currentPeriod = currentMillis - previousMillis;
  while (currentPeriod <= MainPeriod)
  {
    currentMillis = millis();
    currentPeriod = currentMillis - previousMillis;
  }
      
    // need to bufferize to avoid glitches
    unsigned long _duration = duration;
    unsigned long _pulsecount = pulsecount;

    float Freq = 1e6 / float(_duration);  // F = 1/T
    Freq *= _pulsecount; // calculate F

    /*Serial.print(currentMillis);
    Serial.print(" "); // separator!
      //Serial.print(Freq);
    Serial.printf("Freq:%3.2f", Freq);
    Serial.print(" "); 
    Serial.print(_pulsecount);
    Serial.print(" ");
    Serial.println(_duration);
    duration = 0; // clear counters
    pulsecount = 0;*/

    if (Freq > 950.00 && Freq < 1100)
    {
      long sumOfSquares = 0;
      smoothedValue = 0;
      
       for (int i=0; i<numberOfSamples; i++)
       {
          sample = analogRead(A2);         // take a sample reading
          signal = (sample - dc_offset);  // subtract the dc offset to center the signal at 0
          signal *= signal;               // square the value to remove negative values
          sumOfSquares += signal;       // sum the values
        }
  
          // divide the sum of the squares to get the mean
          // and then take the square root of the mean
          RMS = sqrt(sumOfSquares /numberOfSamples); 
  
          // smoothing filter
          //smoothedValue = (alpha * smoothedValue) + ((1-alpha) * RMS);
          
          // convert the RMS value back into a voltage and convert to db  
          //db = 20 * log10(smoothedValue * (aref_voltage / 1024.0)); 
          db = 20 * log10(RMS * (aref_voltage / 1024.0)); 
          dbspl = 94 + db;
          Serial.printf("Freq:%3.2f RMS:%d dB:%-3.2f\n dBSPL:%-3.2f\n", Freq, RMS, db, dbspl);   // output frequency and amplitude data to serial monitor
    }
    else
    Serial.println("Freq detection not in range");
}
    
void myinthandler() // interrupt handler
{
  unsigned long currentMicros = micros();
  duration += currentMicros - previousMicros;
  previousMicros = currentMicros;
  pulsecount++;
}

br
Alan Soong
 
Hi all,

I have make some code for the THD calculation as below.
When I compare the THD% result by using audio analyzer and it is incorrect, can anyone expert help to check whether my code have any issue?

Code: 
void THD1()
{
    //get the current time;
  Starttime = millis();

  while((millis() - Starttime) < Sampletime_ms)
  {
    if (fft1024_1.available())
    {

// each time new FFT data is available
// print it all to the Arduino Serial Monitor
Serial.print("FFT bins ");
for (i = 0; i < 511; i++) {
n = fft1024_1.read(i);

n = n * a_weight[i];
v[i] = n;
v[i] = 0.001953125 *0.3750  * pow(v[i],2); // 1/512 samples = .0019....
magnitude = v[i]+ magnitude;
float p = pow(magnitude,2);

// print only the bins with data for this sine wave binning test
if (n >= 0.01) {
Serial.print(i);
Serial.print(": ");
Serial.print(n);
Serial.print(" ");
Serial.print(p);
Serial.print(" ");
} else {
//Serial.print(" - "); // don't print "0.00"
}
}
Serial.println();

// print the bands (these are NOT ISO bands.....)
//bin size = 44,100 Hz / 1024 or 43.08 Hz per bin
level[0] = fft1024_1.read(0); // FFT "0Hz" -> "DC"
level[1] = fft1024_1.read(1); 
level[2] = fft1024_1.read(23); 
level[3] = fft1024_1.read(46); 
level[4] = fft1024_1.read(69); 
level[5] = fft1024_1.read(93); 
level[6] = fft1024_1.read(116); 

float f1 = 0.001953125 *0.3750 * pow(level[2],2);
float f2 = 0.001953125 *0.3750 * pow(level[3],2);
float f3 = 0.001953125 *0.3750 * pow(level[4],2);
float f4 = 0.001953125 *0.3750 * pow(level[5],2);
float f5 = 0.001953125 *0.3750 * pow(level[6],2);

float p1 = pow(f1,2);
float p2 = pow(f2,2);
float p3 = pow(f3,2);
float p4 = pow(f4,2);
float p5 = pow(f5,2);

Serial.printf("%.6f\n",p1);
Serial.printf("%.6f\n",p2);
Serial.printf("%.6f\n",p3);
Serial.printf("%.6f\n",p4);
Serial.printf("%.6f\n",p5);


        float THD = ((p2 + p3 + p4 + p5)/p1) * 100;
        Serial.printf("THD Percentage:%.6f\n", THD);
        float THD_dB = log10f(THD) * 10.0f;
        Serial.printf("THD_dB:%.2f\n", THD_dB);
     // }
    }
  }
}

Thank you very much!

br,
Alan Soong
 
Alan, it appears that after you find the sum of the powers, magnitude, you apply
Code: 
float p = pow(magnitude,2);
when it should be sqrt.

Also, if you are trying to measure really low values, it may be desirable to omit bins around 44117/24. And also 44117/8. Take a look at the full FFT spectrum to see what is needed. See https://forum.pjrc.com/threads/41356-Audio-Adapter-ADC-Square-wave-with-no-input

Looks like a neat project. Keep us updated.

Bob
 
Hi Bob,

Thank you for your input.

I have added the sqrt to the code and i just take up to 5th harmonic result to calculate the THD.
Unfortunately the result still negative.

From Audio Analyzer, the THD% is about 3%
From my teensy board, the THD% is about 20%.

It is quite big difference.

Thank you and best regards
Alan Soong
 
Hi Alan - I looked further and I believe that when you get to f1, f2, ... you have powers. When you square this with p1, p2,...you have power squared, not what you want.

I suspect the confusion is that in the FFT, Paul takes the square root of the sum of the squares of the in-phase and quadrature outputs, which is the RMS voltage. That is what comes out. Squaring this gives a true estimate of the power in the bin.

An aside is that if projects like yours start running out of real time, going back and redoing the FFT to not take the square root would make things a bit tidier! This is probably not an issue though.

Cheers, Bob

BTW, if you get to working with sine waves not centered on a bin, conservation of energy says you can add the powers of nearby bins to get the true power. Good windowing can make the job easier, as well.
 
Hi Bob,

Is that mean the magnitude = (1/512)* pow(v,2) comes out RMS voltage and it is after square root as indicate in the fft1024.cpp as below?
Code: 
			uint32_t tmp = *((uint32_t *)buffer + i); // real & imag
			uint32_t magsq = multiply_16tx16t_add_16bx16b(tmp, tmp);
			output[i] = sqrt_uint32_approx(magsq);

I have tried to calculate the THD by using f1,f2... instead of p1,p2... and the THD % increase to 70%..

In this case, i have apply Hanning windows and my input source is 1kHz pure tone.

Thanks and best regards
Alan Soong
 
Alan, if I understand correctly, that is not quite right. In your code snippet, magsq is the power in the bin, and output is the rms voltage. Obviously, we are using the terms Voltage and power as if this was an analog circuit---personally this is the easy way to think of these things.

The Hanning sounds great.

Might you provide the whole .ino file and also the output from the Serial Monitor? That might tell a lot.

Also, be sure you are not over-driving the ADC.

Cheers, Bob
 
Hi Bob,

Please see below SPL.ino and also serial output result.PNG for your reference.

Thank you and best regards
Alan Soong
View attachment SPL.ino
serial output result.PNG
 
Alan, in looking at the Serial print, manual calculation of the the second harmonic (adding close bin powers) showed it to be 34% of the total. The software seems OK! Look for an over driving problem. Normally, such issues cause symmetrical clipping that has very little second harmonic, so it is not that simple. Change the level. Check the generator. The readPeakToPeak() of AudioAnalyzePeak object is very useful, but note that a full 16-bits for that is 2.0 output (not 1.0). Good Luck, Bob
 
Hi Bob,

Thanks for your verification.
I have tried on checking the over driving problem and there is no such issue.
In that case, mean my THD software code is correct and applicable?

Thank you and best regards
Alan Soong
 
Audio Function Test source code (With Freq, SPL, THD and SNR)

Hi all,

I'm managed to complete the audio Functional test with functions of freq, SPL, THD and SNR calculation.
I uploaded the source code here so that anyone of you able to use it when needed and also welcome to verify my source code if there is any mistakes.:)

Thank you and best regards
Alan Soong

View attachment AudioFunctionalTest.ino

And My results:
Result.PNG
 
Would it be ok to add this as one of the library's examples? Should it mention your name for credit? A license like MIT or public domain?
 
Hoping the code would show something practical for my use with what was involved ... read it as text. Seemed like a good way to summarize the details and concepts not seen/used before:

SNR (other?) code seems to have a focus hardcoded on reference value 1000 Hz correct?

The comments are nice - but not all accurate? For example::
> level[2] = fft1024_1.read(22,25); // 86- 129Hz -> 108Hz Center
should this actually be::
> level[2] = fft1024_1.read(22,25); // 946- 1075Hz -> 1010Hz Center

SPL works across all bins, but THD and SNR only sample cusps of data? Both skip bin 2-21 then THD skips bins 25-44 (and other groups), and SNR skips bins 26-39. I'm not sure why some bins are excluded from those calculations - it seems they would hold distortion or noise data?

I liked the web link to stackexchangecircuit - was the MCP6022 used? A link to this thread might help find your dev notes.
 
Hi defragster,

Thank you for your comment.
1)For the freq bin comment, I forgot to amend it accordingly. (will update it soon)
2) For THD, I just figured the fundamental till 6th harmonic. Fundamental is around bin 22~ 25, that'w why I skip bin 2~21.
2nd harmonic is around bin 45~48, so in between 26 to 44 is not included on my code.
The rest of the harmonics are also configured with same method, so some group of bins are not included too.
3) For SNR, it is also using almost same method to configure. I'm using power of fundamental signal to divide the Noise power of all harmonic + noise and exclude DC component as well.

I hope that my theory understanding are in correct way. If it is something incorrect, please help to correct me so that I can tidy up these code to make it correct.

Thanks and best regards
Alan Soong
 
Status
Not open for further replies.
Back
Top