Advice on using openmusiclabs fft code with external ADC

Using Arduino Programming Questions
#1

Hello,

I am looking for an FFT library I can use on the Mega2560. Basically, I want to take 256 samples at 391uS intervals. I will end up with a 128 bin spectra from 0-1280 hz I can then discard the last 28 bins and calculate overall vibration acceleration and velocity levels in the 10-1000hz range.

I have an ADC shield I designed with an MCP3208 12 bit ADC. The sensor is an ADXL335 accelerometer giving 300mV/g I am only using the Y axis in the vertical orientation so with no vibration at all, it puts out roughly 1.95V

What I have done so far is to take my samples, find the average, which as I mentioned is about 1.95V which I subtract ie remove the dc bias, convert them to g’s, pass them through a hanning window then do an fft with some old code I adapted from an old Turbo C++ program I developed years ago.

Now, as the project is developing and I’m adding in more bits, I find that the old fft routine is causing problems.

OK, the openmusiclabs fft code takes its data from the internal ADC and looks like this

      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

and my code, using the MCP3208 library, looks like this

void sample(int chan) 
{
  unsigned long start_time;
  unsigned long next_time;
  
  int a;
  start_time = micros();
  if(start_time>0xFFF0BDBF) // micros is close to rollover, wait until it has.
    delay(1010);            // 1.01 second delay to take micros past rollover

  next_time = micros()+391;
  for(a=0; a<=numsamples; a++)
  {
    twf[a] = adc.analogRead(chan)*res; // res=3.3/4096 ie 12 bit resolution with 3.3V ref
    avg+=twf[a];

    for(;;)
      if((long)micros()>next_time)
        break;

  next_time+=391;
  }
  avg/=numsamples;
}

So, I’m guessing that ADCL and ADCH are the low and high bytes as the internal ADC is 10 bits and they are shifted into int k and fft_input is made up of alternative real and imaginary values but I’d like to be corrected if I’m wrong. Also, this code does not remove any dc bias prior to doing the fft which I have seen on other examples.

Not being a brilliant mathematician (or a brilliant anything for that matter) I would really appreciate some advice here. Is it necessary to convert the voltages to their relative values prior to doing the fft? should the dc offset be removed? or can I just fetch the data from the MCP3208 and put it into the even bins, fill the odd bins with zeros, do the fft and worry about conversion later?

I know there’s some knowledgeable people on here who seem to understand fft’s very well as I’ve been reading your posts all afternoon and I’d be very grateful for your advice.

#2 
      k -= 0x0200; // form into a signed int
      k <<= 6; // form into a 16b signed int

The comments do not match the action

      k -= 0x0200; // subtract offset (512)
      k <<= 6; // multiply result by 64
#3

CM_Guy:
So, I'm guessing that ADCL and ADCH are the low and high bytes ...

Indeed. No need to guess, though. You can find that information in the datasheet, linked on the Uno's product page, here - Arduino Uno Rev3 — Arduino Official Store. Chapter 24 tells you all about the ADC, including the registers associated with it.

... the internal ADC is 10 bits and they are shifted into int k and fft_input is made up of alternative real and imaginary values...

Yes.

Also, this code does not remove any dc bias prior to doing the fft which I have seen on other examples.

~~No, it doesn't.~~It does, if the DC offset is at midscale on a ten-bit ADC. With your DC offset at some other value, and a 12-bit ADC, you'll need to do something different. [Edit: corrected this response, in accordance with el_supremo's observation, below.]

Is it necessary to convert the voltages to their relative values prior to doing the fft?

No. You can make linear conversions before you run it, or you can make them after you run it. You can verify that by looking up the definition of the discrete Fourier transform, and consider what happens when you multiply each data term by a constant.

should the dc offset be removed?

You certainly can. If you're using a Hanning window, removing the DC offset will likely improve your ability to resolve peaks in low-numbered bins.

#4

tmd3:
Indeed. No need to guess, though.

Hello tmd3. Thank you for your help. It is nice to know that I'm thinking along the right lines. Now, the reason I was guessing is not because I cant be arsed reading chapter24 to find out about the ADC, but I don't see the need. I put my time into learning how to use the MCP3208 as 12 bits are the minimum I need for my work. Its like me learning to speak Inuit when I really need Greek because I live in UK and Cyprus! However, thank you for the link because it may well point other people in the right direction where 10 bits IS sufficient.

I think that you have given me enough information for me to give it a go and see what output I get.

@Whandall, thanks for the post. However, that is how the code reads. I just pasted it in direct from their code but I agree, its confusing! btw, I prefer 'so long and thanks for all the fish' but then, I just love dolphins. We occasionally get pods of them round the yacht club.

#5

CM_Guy:
I just pasted it in direct from their code but I agree, its confusing!

The comments are plain wrong, not confusing.

#6 
      int k = (j << 8) | m; // form into an int
      k -= 0x0200; // form into a signed int
      k <<= 6; // form into a 16b signed int

This code does do the DC offset. The first statement forms the high and low bytes into a signed integer (as long as int means 16-bit integer - if int is 32-bits this won’t work).
The second statement subtracts 512 which changes the ADC reading from a range of 0-1023 to -512 to 511. This is the DC offset.
The third statement shifts the number up 6 bits which, in effect, amplifies the signal by a factor of 64. This will also have the effect of amplifying the noise and adding more quantization noise. I’m not sure that it is useful to do this but the rest of the code might rely on it.

Pete

#7

el_supremo:
This code does do the DC offset.

Indeed it does, for a specific offset. That post has been corrected. It's useful if the actual DC offset is half-scale on a ten-bit ADC, or otherwise corresponds to 512 ADC counts. Given that the traditional offset scheme for an Arduino is a voltage divider at midscale, with the AC signal capacitively coupled, that's certainly a reasonable default for example code.

A more general DC offset would take the average value of the signal over all the samples, and subtract that from each sample. An alternative might be to subtract the values of the windowed samples, to drive the DC term from the FFT to zero. Windowing tends to spread the DC term over the low-numbered bins - a low-amplitude, low-frequency signal can be swamped by spectral leakage from bin zero. As to whether it's more theoretically sound to eliminate DC before or after windowing, that's beyond me.

The third statement shifts the number up 6 bits which, in effect, amplifies the signal by a factor of 64.

I believe that the sample values are adjusted in order to reduce the effect of rounding in the intermediate calculations of the FFT. For space and speed, those are done with fixed-point math, and each result is adjusted to a 16-bit integer, so a little information is rounded away in each step. Starting with a larger signal value reduces the effect of that rounding.

#8

Many thanks for your input guys. I do have it working quite nicely. So the sampling code looks like this

//Function fftsample(int chan, int b)
// reads samples from channel chan into global array fft_input
// b is step. use b=1 to get twf data for export
// use b=2 to get alternative twf data and zero's for fft
// ie real and imaginary data are interspersed

void fftsample(int chan, int samples, int b)
{
     unsigned long start_time;
     unsigned long next_time;
     double            twftot=0.0;
     int                  a;
  
    if(micros()>0xFFF0BDBF)   // micros is less than 1 second until rollover
      delay(1010);            // 1.01 second delay to take micros past rollover
    
    start_time = micros();
    next_time = micros()+391;
    for(a=0; a<samples; a+=b)
    {
        fft_input[a] = adc.analogRead(chan);    // put real data into even bins
        twftot+=fft_input[a];
        fft_input[a+1] = 0;                             // set odd bins to 0
    
        for(;;)
            if((long)micros()>next_time)
                break;

        next_time+=391;
    }

    twfavg=(int) (twftot/(samples/b));

}

and in the loop, the fft bit is

fftsample(p,512,2);
            
for(i=0; i<512; i+=2)
    fft_input[i]-=twfavg;

          
                       
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
            
for(i=4; i<200; i+=2)
{
    fft_input[i]=absolute(fft_input[i],fft_input[i+1]);
    acceleration=fft_input[i]*res*(1000/sensitivity);
    acum+=(acceleration*acceleration);
    velocity=(acceleration*9806.65)/(twopi*10*i);
    vcum+=(velocity*velocity);
}

acc[p]=sqrt(acum/lineshape);
vel[p]=sqrt(vcum/lineshape);

So in the sampling part, I calculate the average which is then subtracted from the twf prior to ding the fft. I then calculate the power in both the acceleration and velocity spectrums and go on to display the results.

Next question, has anyone modified the openmusiclabs fft code to allow larger fft’s - eg 512 samples rather than the maximum 256 as provided. That would just be the icing on today’s cake.

Thanks again to all. Steve.

#9

The openmusiclabs FFT code is intended to run on an Uno. The Uno has 2K SRAM, and, with integer data at two bytes per sample, a 256-point array, together with the imaginary components, takes 1K. The FFT code requires a number of samples that's an integer power of two, so, the next increment of increase would take 2K, or all of the SRAM, leaving nothing for the incidental pointers and counters that control the calculation. That's why the library doesn't support FFT's bigger than 256 samples.

The FFT library uses precaluclated twiddle factors, stored in program memory. The library provides those for no more than 256 samples. It would be generally easy to establish new twiddle factors for 512 samples, but I don't know that it would help. Note that the heart of the FFT library is written in assembler, and doesn't easily yield to analysis or modification.

The openmusiclabs FFT library is written for speed. If your application will allow for a slower calculation, you can certainly implement an FFT function of your own, and its limits can be whatever you like, and can manage with the 8K of SRAM in the Mega 2560. There are a whole bunch of implementations available on the web, written in C or C++. Some are already crafted for the Arduino, and some are more general, and will require some effort to massage them into your code. It may be hard to beat the speed of the openmusiclabs code on an Arduino, but it's not hard to beat the size.

#10

As to whether it's more theoretically sound to eliminate DC before or after windowing, that's beyond me.

It's quite simple. A DC offset before, will limit the dynamic range of the calculations by lowering the clip level.

#11