Forum Rule: Always post complete source code & details to reproduce any issue!
Results 1 to 17 of 17

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

  1. #1

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

    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

  2. #2
    Senior Member PaulStoffregen's Avatar
    Join Date
    Nov 2012
    Posts
    22,079
    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?

  3. #3
    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

  4. #4
    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

  5. #5
    Senior Member
    Join Date
    Dec 2016
    Posts
    118
    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...-with-no-input

    Looks like a neat project. Keep us updated.

    Bob

  6. #6
    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

  7. #7
    Senior Member
    Join Date
    Dec 2016
    Posts
    118
    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.

  8. #8
    Hi Bob,

    Is that mean the magnitude = (1/512)* pow(v[i],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

  9. #9
    Senior Member
    Join Date
    Dec 2016
    Posts
    118
    Alan, if I understand correctly, that is not quite right. In your code snippet, magsq is the power in the bin, and output[i] 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

  10. #10
    Hi Bob,

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

    Thank you and best regards
    Alan Soong
    SPL.ino
    Click image for larger version. 

Name:	serial output result.PNG 
Views:	43 
Size:	16.1 KB 
ID:	9712

  11. #11
    Senior Member
    Join Date
    Dec 2016
    Posts
    118
    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

  12. #12
    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

  13. #13

    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

    AudioFunctionalTest.ino

    And My results:
    Click image for larger version. 

Name:	Result.PNG 
Views:	75 
Size:	15.6 KB 
ID:	9871

  14. #14
    Senior Member PaulStoffregen's Avatar
    Join Date
    Nov 2012
    Posts
    22,079
    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?

  15. #15
    Hi Paul,

    It is fine for me.
    It can used for public domain.

    Thank you and best regards
    Alan Soong

  16. #16
    Senior Member+ defragster's Avatar
    Join Date
    Feb 2015
    Posts
    11,764
    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.

  17. #17
    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

Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts
  •