Incorrect/inconsistent frequency readings for clarinet sound

Hello everyone!
I am working on a project for my music professor that is sort of a musical sculpture: a derelict clarinet has different colored LEDs in its tone holes and each one lights depending on the frequency of the note played by a separate instrument.

Of course, for this I need to detect the frequency of the played note with reasonable accuracy, and the problem I am having is that with one method (Amanda Ghassei’s zero crossing-esque algorithm), I get wildly inconsistent readings (to the point where the clarinet’s high C, which should read at ~1000 Hz is reading at 100), and with the other (OpenMusicLabs FFT library), I a) don’t seem to know enough to correctly interpret the output and b) that is not helped by the fact that almost every “bin” from the FFT output is similar in value.

Therefore, I need help with a) understanding the FFT library enough to use it (preferred, as this appears to be the “correct” way) or b) figuring out where I’m going wrong with my implementation of Amanda Ghassei’s algorithm.

My setup is as follows:
I am detecting sound by means of an electret microphone connected to a TLC272 op-amp and biased around 2.5V (schematic attached). That output is connected to Analog 0, where it is read by the code.

Code will be posted in the first reply to get around the post-length limit.

Thank you for your help!

Here is my adaptation of the zero-crossing-esque algorithm:

//Experimentation code for frequency detection 
//Asa Graf 11/27/2014



//Adapted from: 

//generalized wave freq detection with 38.5kHz sampling rate and interrupts
//by Amanda Ghassaei
//http://www.instructables.com/id/Arduino-Frequency-Detection/
//Sept 2012

/*
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or
 * (at your option) any later version.
 *
*/

int indicatorPin = 13;
int brightPin = 3;

#define filterSamples 10
float frequencyAverage = 0;
  
int readings[filterSamples];


//clipping indicator variables
boolean clipping = 0;

//data storage variables
byte newData = 0;
byte prevData = 0;
unsigned int time = 0;//keeps time and sends vales to store in timer[] occasionally
int timer[10];//sstorage for timing of events
int slope[10];//storage for slope of events
unsigned int totalTimer;//used to calculate period
unsigned int period;//storage for period of wave
byte index = 0;//current storage index
float frequency;//storage for frequency calculations
int maxSlope = 0;//used to calculate max slope as trigger point
int newSlope;//storage for incoming slope data

//variables for decided whether you have a match
byte noMatch = 0;//counts how many non-matches you've received to reset variables if it's been too long
byte slopeTol = 3;//slope tolerance- adjust this if you need
int timerTol = 10;//timer tolerance- adjust this if you need

//variables for amp detection
unsigned int ampTimer = 0;
byte maxAmp = 0;
byte checkMaxAmp;
byte ampThreshold = 30;//raise if you have a very noisy signal

void setup(){
  
  Serial.begin(9600);
  
  pinMode(indicatorPin,OUTPUT);//led indicator pin
  pinMode(12,OUTPUT);//output pin
  
  pinMode(brightPin, OUTPUT); //LED brightness pin
  
  cli();//diable interrupts
  
  //set up continuous sampling of analog pin 0 at 38.5kHz
 
  //clear ADCSRA and ADCSRB registers
  ADCSRA = 0;
  ADCSRB = 0;
  
  ADMUX |= (1 << REFS0); //set reference voltage
  ADMUX |= (1 << ADLAR); //left align the ADC value- so we can read highest 8 bits from ADCH register only
  
  ADCSRA |= (1 << ADPS2) | (1 << ADPS0); //set ADC clock with 32 prescaler- 16mHz/32=500kHz
  ADCSRA |= (1 << ADATE); //enabble auto trigger
  ADCSRA |= (1 << ADIE); //enable interrupts when measurement complete
  ADCSRA |= (1 << ADEN); //enable ADC
  ADCSRA |= (1 << ADSC); //start ADC measurements
  
  sei();//enable interrupts
}

ISR(ADC_vect) {//when new ADC value ready
  
  PORTB &= B11101111;//set pin 12 low
  prevData = newData;//store previous value
  newData = ADCH;//get value from A0
  if (prevData < 127 && newData >=127){//if increasing and crossing midpoint
    newSlope = newData - prevData;//calculate slope
    if (abs(newSlope-maxSlope)<slopeTol){//if slopes are ==
      //record new data and reset time
      slope[index] = newSlope;
      timer[index] = time;
      time = 0;
      if (index == 0){//new max slope just reset
        PORTB |= B00010000;//set pin 12 high
        noMatch = 0;
        index++;//increment index
      }
      else if (abs(timer[0]-timer[index])<timerTol && abs(slope[0]-newSlope)<slopeTol){//if timer duration and slopes match
        //sum timer values
        totalTimer = 0;
        for (byte i=0;i<index;i++){
          totalTimer+=timer[i];
        }
        period = totalTimer;//set period
        //reset new zero index values to compare with
        timer[0] = timer[index];
        slope[0] = slope[index];
        index = 1;//set index to 1
        PORTB |= B00010000;//set pin 12 high
        noMatch = 0;
      }
      else{//crossing midpoint but not match
        index++;//increment index
        if (index > 9){
          reset();
        }
      }
    }
    else if (newSlope>maxSlope){//if new slope is much larger than max slope
      maxSlope = newSlope;
      time = 0;//reset clock
      noMatch = 0;
      index = 0;//reset index
    }
    else{//slope not steep enough
      noMatch++;//increment no match counter
      if (noMatch>9){
        reset();
      }
    }
  }
    
  if (newData == 0 || newData == 1023){//if clipping
    PORTB |= B00100000;//set pin 13 high- turn on clipping indicator led
    clipping = 1;//currently clipping
  }
  
  time++;//increment timer at rate of 38.5kHz
  
  ampTimer++;//increment amplitude timer
  if (abs(127-ADCH)>maxAmp){
    maxAmp = abs(127-ADCH);
  }
  if (ampTimer==1000){
    ampTimer = 0;
    checkMaxAmp = maxAmp;
    maxAmp = 0;
  }
  
}

void reset(){//clea out some variables
  index = 0;//reset index
  noMatch = 0;//reset match couner
  maxSlope = 0;//reset slope
}


void checkClipping(){//manage clipping indicator LED
  if (clipping){//if currently clipping
    PORTB &= B11011111;//turn off clipping indicator led
    clipping = 0;
  }
}


void loop(){
  
  checkClipping();
  //Average 10 readings of frequency 
  
  if (checkMaxAmp>ampThreshold){
    frequency = 38462/float(period);//calculate frequency timer rate/period
    int f = (int) frequency;
    
    int fSmooth = digitalSmooth(f, readings);
    
    Serial.print(frequency);Serial.print("               ");//Serial.print(fSmooth);
    //Serial.println(" hz");
    int note = map(frequency, 100,1000,1,33);
    Serial.println(note);
    int brightness = constrain(map(note, 1, 33, 0,255),0,255);

    analogWrite(brightPin,brightness);

    delay(1);
  }
  

}




int digitalSmooth(int rawIn, int *sensSmoothArray){     // "int *sensSmoothArray" passes an array to the function - the asterisk indicates the array name is a pointer
  int j, k, temp, top, bottom;
  long total;
  static int i;
 // static int raw[filterSamples];
  static int sorted[filterSamples];
  boolean done;

  i = (i + 1) % filterSamples;    // increment counter and roll over if necc. -  % (modulo operator) rolls over variable
  sensSmoothArray[i] = rawIn;                 // input new data into the oldest slot

  // Serial.print("raw = ");

  for (j=0; j<filterSamples; j++){     // transfer data array into anther array for sorting and averaging
    sorted[j] = sensSmoothArray[j];
  }

  done = 0;                // flag to know when we're done sorting              
  while(done != 1){        // simple swap sort, sorts numbers from lowest to highest
    done = 1;
    for (j = 0; j < (filterSamples - 1); j++){
      if (sorted[j] > sorted[j + 1]){     // numbers are out of order - swap
        temp = sorted[j + 1];
        sorted [j+1] =  sorted[j] ;
        sorted [j] = temp;
        done = 0;
      }
    }
  }

/*
  for (j = 0; j < (filterSamples); j++){    // print the array to debug
    Serial.print(sorted[j]); 
    Serial.print("   "); 
  }
  Serial.println();
*/

  // throw out top and bottom 15% of samples - limit to throw out at least one from top and bottom
  bottom = max(((filterSamples * 15)  / 100), 1); 
  top = min((((filterSamples * 85) / 100) + 1  ), (filterSamples - 1));   // the + 1 is to make up for asymmetry caused by integer rounding
  k = 0;
  total = 0;
  for ( j = bottom; j< top; j++){
    total += sorted[j];  // total remaining indices
    k++; 
    // Serial.print(sorted[j]); 
    // Serial.print("   "); 
  }

//  Serial.println();
//  Serial.print("average = ");
//  Serial.println(total/k);
  return total / k;    // divide by number of samples
}

As you can see, I attempted a sort of digital smoothing algorithm but, while that helped to smooth out outlier measurements, it didn’t change the fact that sometimes I read completely incorrect frequencies.

Here is the code that I have been using for the FFT: It is merely one of the examples included in the library.

/*
fft_adc_serial.pde
guest openmusiclabs.com 7.7.14
example sketch for testing the fft library.
it takes in data on ADC0 (Analog0) and processes them
with the fft. the data is sent out over the serial
port at 115.2kb.
*/

#define LOG_OUT 1 // use the log output function
#define FFT_N 256 // set to 256 point fft

#include <FFT.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 < 512 ; i += 2) { // 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
      fft_input[i] = k; // put real data into even bins
      fft_input[i+1] = 0; // set odd bins to 0
    }
    fft_window(); // window the data for better frequency response
    fft_reorder(); // reorder the data before doing the fft
    fft_run(); // process the data in the fft
    fft_mag_log(); // take the output of the fft
    sei();
    Serial.println("start");
    for (byte i = 0 ; i < FFT_N/2 ; i++) { 
      Serial.println(fft_log_out[i]); // send out the data
    }
  }
}

(schematic attached).

You could have fooled me. Such a huge picture ( have you clicked on it? ) with so little in it, pictures should be no more than 1000 pixels wide.

What you are trying to do is almost impossible. This is because the harmonic contents of the note change over time and it is not always the fundamental that is the strongest harmonic. Hence inconstant readings are to be expected. You can get the note but the octave is going to be hit and miss. With a much better processor than the Arduino it is possible to get close but it will still not be perfect.

I do not think zero crossing pitch detection is adequate for more complex waveforms. Of course it will work for all schoolbook exapmles of sine, square sawtoth and so on, but i have never seen these waveforms in real life. So i suggest a fft based method. A arduino uno has very little processing power but it can be done, have a look here. A arduino due or a teensy 3 is much better suited to this task though.

As you can see, I attempted a sort of digital smoothing algorithm but, while that helped to smooth out outlier measurements, it didn’t change the fact that sometimes I read completely incorrect frequencies.

You may need to throw-out the outliers. You know some frequencies (like 100Hz) are entirely invalid if they are beyond the instruments range or in-between notes… You should know what (fundamental) frequencies/ the clarinet can play and what frequencies are noise, harmonics, and overtones, or just erroneous measurements/calculations… and your software should “know” that too!

It would be a good idea if you can get hold of a “real” spectrum analyzer so you know what you’re dealing with. If you are at a university, there should be a frequency analyzer available in the electronics or physics departments if the music department doesn’t have one. Or, you can find software frequency analyzer.

A clarinet shouldn’t be too difficult (compared to other instruments) because it can’t play chords and its not dynamic (the loudness doesn’t vary much). But some instruments (such as lowest notes on the piano), have will have harmonics that are stronger than the fundamental.* Your ear will pick-out the fundamental, but it’s not the dominate pitch/frequency.

* P.S.
I found [u]this paper[/u] that shows some clarinet harmonics are also stronger (or nearly as strong) than the fundamental. I’m not sure if that paper shows the frequency analysis of all of the clarinet notes, but if you can find that information you won’t need to do it yourself with a spectrum analyzer.

I agree that the fft should work much better than a zero-crossing algorithm.

The fft_mag_log() function in the OpenMusicLabs library is intended to give a representation of volume levels as we hear them. It seriously distorts the amplitude scale, so use fft_mag_lin() instead.