Adafruit_ZeroFFT on Feather RP2040--bug?

Please tell us which board you are using.
For CircuitPython issues, ask in the Adafruit CircuitPython forum.

Moderators: adafruit_support_bill, adafruit

Please be positive and constructive with your questions and comments.
Locked
User avatar
cf20855
 
Posts: 9
Joined: Sun Sep 26, 2010 11:07 pm

Adafruit_ZeroFFT on Feather RP2040--bug?

Post by cf20855 »

I'm trying to make a little vibration analyzer using an ADXL345 connected to a Feather RP2040 running FFT software. The idea is to use the ADXL345 to sense vibration and the RP2040 to compute the amplitude of the vibration (Y-axis) vs frequency (X-axis). For test purposes, to generate a single-frequency vibration source, I drive my stereo system with a sine wave output by a signal generator. I place the ADXL345 on the woofer speaker cone. I pipe the FFT outputs out of the Feather RP2040 to a Python program on my laptop to graph the results on an amplitude vs. frequency plot.

I first collected the vibration samples and performed the FFT computation using a CircuitPython program on the Feather RP2040. I used spectrogram from the ulab.scipy.signal library to compute the amplitude of the vibration data vs. frequency. That seemed to work pretty well, but I wanted a speed improvement. I decided to give the Adafruit_Zero_FFT_Library a try, under the Arduino IDE 2.0 environment.

When I plotted the FFTZero output, however, I noticed a strange phenomenon. The plotted signal amplitude slowly wandered up and down. This had not happened with the CircuitPython program. I modified my data collection routine to print the measured acceleration amplitude samples instead of the amplitude vs frequency output of the FFT, and the FFT input looked fine.

I decided to delve into the fftutil.c code in the Adafruit_Zero_FFT_Library, and noticed something that seemed odd to me at the bottom of the code:

Code: Select all

for (int i = 0; i < length; i++) {
	q15_t val = *pOut++;
	uint32_t v = abs(val);
	pSrc++ = v;
	pOut++; // discard imaginary phase val
}
I have a very lightweight background in digital signal processing, but my understanding is that if one wants an FFT algorithm to output amplitude vs. frequency, one takes the square root of the sum of the squares of the real and complex outputs, and doesn't just throw away the complex vector. This is what spectrogram from the ulab.scipy.signal library does--see https://micropython-ulab.readthedocs.io ... pectrogram. Moreover, the complex values aren't the phases of the signal in the frequency bins. Rather, the phase of the signal in a particular frequency bin is the arctangent of the complex value over the real value in that bin. And the phase of the ADXL345 sampling clock varies with respect to the phase of the sine wave driving my woofer, which I thought could explain the results I saw from this code.

As an aside, the terms "real" and "complex" vectors seem a bit inapt. Perhaps a better understanding the FFT output vectors would be obtained by calling them the "real" and "imaginary" components of the signals in the frequency bins. Refer to https://www.dsprelated.com/showarticle/192.php for example.

As another aside, the square root computation is a little bit expensive, but there are some very good magnitude approximations that do not require taking the square root of the sum of the squares of the real and complex values. Such approximations are documented at https://dspguru.com/dsp/tricks/magnitude-estimator/, for example, and in a paper titled "A Baker's Dozen Magnitude Approximations and Their Detection Statistics".

I modified the above code to read as shown below, using one of the magnitude approximation formulas, adjusted to q15 arithmetic. It seems to work OK at first glance, but more testing is required.

Code: Select all

for (int i = 0; i < length; i++) {
	// r denotes real, i denotes imaginary
	q15_t r_val = *pOut++;
	q15_t i_val = *pOut++;
	uint32_t r_v = abs(r_val);
	uint32_t i_v = abs(i_val);
	// Mag ~=Alpha * max(|r|, |i|) + Beta * min(|r|, |i|)
	// where Alpha = 0.960433870103 and Beta = 0.397824734759 (in q15 representation)
	uint32_t v = (31471 * max(r_v, i_v) + 13036 * min(r_v, i_v)) >> 15;
	*pSrc++ = v;
}

User avatar
adafruit_support_mike
 
Posts: 67391
Joined: Thu Feb 11, 2010 2:51 pm

Re: Adafruit_ZeroFFT on Feather RP2040--bug?

Post by adafruit_support_mike »

Your information about quadrature signals is valid, but mostly irrelevant to FFT calculations.

Quadrature eliminates Nyquist Folding, where a 999Hz signal and a 1001Hz signal sampled at 1000Hz both alias to 1Hz.. strictly speaking the 999Hz signal aliases to -1Hz while the 1001Hz signal aliases to 1Hz, but there's no way to tell those frequencies apart just looking at the waveforms. Quadrature makes it possible to distinguish positive frequencies from negative frequencies.

(Aliasing to a negative frequency is what makes wagon wheels in old movies look like they're rotating slowly backwards)

FFTs should never see a collection of samples where negative frequencies are possible. That breaks the Nyquist criterion that the sampling frequency must be more than twice the highest frequency of interest to reconstruct frequency and amplitude. You can identify the existence of most signals at exactly half the sampling frequency, but the sampler will always catch the waveform at the same level, making the signal alias to a DC value that depends on the relative phases of the signal and the sample collector.

To get frequency and amplitude, you need more than two samples per cycle of the input waveform so the sample collector will (eventually) see the input signal at all amplitudes.

Now, it is common to represent periodic functions with complex exponentials.. Euler's great insight was that doing so reduces trigonometry to algebra and makes life a heck of a lot easier. When we look at complex exponentials as a mental model, they represent a point rotating around the X axis, at frequency equal to the x coordinate, in an X-Y-i coordinate space.

FFTs don't care about that, though. They operate on sine waves.. specifically the fact that the product of sine waves at two different frequencies is equivalent to the sum of sines at two other frequencies: (f1+f2)+(f1-f2).

That waveform, being composed of two sine waves, spends equal amounts of time positive and negative of the X axis, so its integral over a sufficiently long time will be zero.

There's one catch though: if you multiply a sine wave by itself, every output value is the square of some input value. Squares can't be negative, so the integral of a sine wave multiplied by itself is proportional to the amplitude of the input waveforms.

The Fourier transform does that 'multiply by frequency f(x) then take the integral' with an arbitrary waveform as one input and a sine wave at every frequency in a given range as the other input. For any value of x, the product of the input and f(x) turns any amount of f(x) in the input to a value that can never be negative, and components of every other frequency into values that will integrate to zero. Therefore any nonzero value in the integral must come from an f(x) component of the input.

The graph of all those integrals tells you all the frequency components that make up the input.

The problem with the Fourier transform is that it's a continuous function and extremely difficult to calculate mechanically.

The FFT solves a simpler version of the problem: discrete Fourier transforms. It uses a series of samples as one input, and values of f(x) sampled at the same frequency, for some discrete set of f(x) frequencies. That makes all the calculations straightforward, but the number of calculations you have to do increases with the square of the number of samples. Programmers refer to that as O(n^2) time, which gets prohibitively expensive no matter how fast your hardware is.. eventually we can find a value of n whose n^2 is too big to be reasonable.

One property of the regular DFT is that you end up having to re-calculate products of the same values over and over. In 1965 James Cooley and John Tukey worked out a way to calculate products once, save them, and simply look them up when they're needed again. That reduced the amount of work to O(n*log(n)), which is low enough to be reasonable.. now almost every complicated signal in the world gets FFT'd somewhere along the way.


Getting back to complex numbers and frequencies: the complex/imaginary part of an exponential frequency shows how signals would behave if they were helices in a complex space, but that's an artifact of the tool we used to make the math easier, not a physical reality. When we want to talk about the actual sine wave, we drop the complex/imaginary component and only pay attention to the real component.

cf20855 wrote: Mon Sep 19, 2022 10:50 am When I plotted the FFTZero output, however, I noticed a strange phenomenon. The plotted signal amplitude slowly wandered up and down.
What value are you referring to? An FFT delivers a table of amplitudes that represent frequency components of the input.

If you're just selecting the amplitude at a specific frequency, the up-and-down suggests your selected frequency is slightly off from the frequency in the samples.. you're capturing part of the (f1+f2)+(f1-f2) where f1 and f2 are close together, but not exactly equal.

One side effect of the DFT is that it can't identify frequencies exactly, since it only multiplies-and-integrates a discrete set of reference frequencies. Real frequency components near each other get assigned to 'bins' in the output.

Try plotting the bins on either side of the one you're using now, and see if they show an equal-but-opposite wobble.

User avatar
cf20855
 
Posts: 9
Joined: Sun Sep 26, 2010 11:07 pm

Re: Adafruit_ZeroFFT on Feather RP2040--bug?

Post by cf20855 »

To clarify my experiment: I generated single-frequency tones using a signal generator, amplified them, and used a woofer to translate the electrical signals to mechanical vibration. I varied the tone frequencies over a 10-200 Hz range, and could either see the sub-20 Hz motion of the woofer or hear the sound it produced at higher frequencies. I used the ADXL345 to sample the mechanical vibration, holding it on the speaker cone. I set the ADXL345 sample rate to 400 samples/second, yielding a 200 Hz bandwidth. I collected groups of 256 vibration samples as the input to the FFT algorithm, resulting in amplitude values in 128 frequency bins. I plotted the amplitude vs frequency output of the FFT at 128 frequencies from 0-200 Hz.

I saw the amplitude variation occurring on single-frequency sinusoidal signals across at least a 30-125 Hz frequency range, when using the unmodified Adafruit library code. I did not see this amplitude variation when using my modified code. The CircuitPython code, which uses a floating point FFT, produced stable amplitude vs frequency data over at least a 10 Hz to 195 Hz range.

You can easily demonstrate the issue with the FFTZero library I pointed out in my original post using a modified version of the fft_test example included in the Adafruit ZeroFFT library:

Code: Select all

/* This example shows the most basic usage of the Adafruit ZeroFFT library.
 * it calculates the FFT and prints out the results along with their corresponding frequency
 * 
 * The signal.h file constains a 200hz sine wave mixed with a weaker 800hz sine wave.
 * The signal was generated at a sample rate of 8000hz.  This vector is called 'asignal'.
 * The vector 'bsignal' is a randomly phase-shifted sample of 'asignal', with a length of DATA_SIZE.
 * 
 * Note that you can print only the value (comment out the other two print statements) and use
 * the serial plotter tool to see a graph.  This version of fft_test.ino is modified to
 * randomly vary the sampling phase with respect to the signal phase, do the FFT, and 
 * output the FFT output values in a loop.
 */

#include "Adafruit_ZeroFFT.h"
#include "asignal.h"

//the signal in signal.h has 2048 samples. Set this to a value between 16 and 2048 inclusive.
//this must be a power of 2
#define DATA_SIZE 256

q15_t bsignal[DATA_SIZE];

//the sample rate
#define FS 8000

// the setup routine runs once when you press reset:
void setup() {
  Serial.begin(115200);
  while(!Serial); //wait for serial to be ready

  //run the FFT
  // ZeroFFT(asignal, DATA_SIZE);

  //data is only meaningful up to sample rate/2, discard the other half
  // for(int i=0; i<DATA_SIZE/2; i++){
    
    //print the frequency
    //Serial.print(FFT_BIN(i, FS, DATA_SIZE));
    //Serial.print(" Hz: ");

    //print the corresponding FFT output
    // Serial.println(asignal[i]);
  //}
}

void loop() {
  while(true) {

    int offset = int(random(DATA_SIZE));

    for(int i=0; i<DATA_SIZE; i++){
      bsignal[i] = asignal[(i + offset) % 2048];
    }

    //run the FFT
    ZeroFFT(bsignal, DATA_SIZE);

    //data is only meaningful up to sample rate/2, discard the other half
    for(int i=0; i<DATA_SIZE/2; i++){
      
      //print the frequency
      //Serial.print(FFT_BIN(i, FS, DATA_SIZE));
      //Serial.print(" Hz: ");

      //print the corresponding FFT output
      Serial.print(bsignal[i]);
      Serial.print(",");
    }
  Serial.println();
  delay(1000);
  }
}

Try it with the original version of fftutil.c and with a version of fftutil.c modified as I show in my original post. You can plot the output of fft_test.ino using a little Python program I have tailored for this purpose:

Code: Select all

# This program was originally written to plot the output of an FFT performed on accelerometer data
import serial
import matplotlib.pyplot as plt
from drawnow import *
import atexit
import numpy as np

values = []
valuesReadArray = []

plt.ion()
cnt=0
values = np.zeros(128)

# Modify the com port setting in the following line of code to match what the Arduino IDE is using
serialAccel = serial.Serial('COM3', 115200)
print("serialAccel.isOpen() = " + str(serialAccel.isOpen()))


def plotValues():
    ax = plt.gca()
    plt.title('FFTZero Output for Varying Phases of asignal')
    plt.grid(True)
    plt.xlabel('Hertz')
    plt.ylabel('Values')
    
    # Set y axis range.
    ax.set_ylim(-10, 2000)
    
    # Set x axis tick labels
    ax.set_xticklabels(["0", "0", "625", "1250", "1875", "2500", "3125", "3750"])

    plt.plot(values, 'rx-', label='values')
    plt.legend(loc='upper right')

def doAtExit():
    serialAccel.close()
    print("Close serial")
    print("serialAccel.isOpen() = " + str(serialAccel.isOpen()))

atexit.register(doAtExit)

while True:
    while (serialAccel.inWaiting()==0):
        pass

    valuesRead = serialAccel.readline(1500)
    
    valuesReadStr = str(valuesRead)

    valuesReadArray = valuesRead.split(b',')
    
    lengthArray = len(valuesReadArray)
    
    for i in range(1, lengthArray-2):
        values[i] = int(valuesReadArray[i])
        
    drawnow(plotValues)
And in case anyone is interested, here is the CircuitPython code that samples the accelerometer data and performs an FFT on it:

Code: Select all

# SPDX-FileCopyrightText: 2021 ladyada for Adafruit Industries
# SPDX-License-Identifier: MIT
# Vibration spectrogram by Chip F. 11 September 2022

import adafruitcf_adxl34x
import board
import busio
import time


from ulab import numpy as np
from ulab.scipy.signal import spectrogram

# The following I2C initialization results in a SCL frequency of about 400 kHz... YMMV
i2c = busio.I2C(board.GP3, board.GP2, frequency=500_000)

# For ADXL345
accelerometer = adafruitcf_adxl34x.ADXL345(i2c)

# Configure the ADXL345 sampling rate and g range
accelerometer.data_rate = adafruitcf_adxl34x.DataRate.RATE_400_HZ
accelerometer.range = adafruitcf_adxl34x.Range.RANGE_8_G

# Set the FFT length
fft_size = 256

# Initialize arrays for x, y, and z samples and the Blackman window array
samples = np.zeros((fft_size, 3))
samples_x = np.zeros(fft_size)
samples_y = np.zeros(fft_size)
samples_z = np.zeros(fft_size)
bw = np.zeros(fft_size)

# Compute the Blackman Window coefficients
for i in range(fft_size):
    bw[i] = (
        0.42
        - 0.5 * np.cos(2 * np.pi * i / fft_size)
        + 0.08 * np.cos(4 * np.pi * i / fft_size)
    )
    
# Loop forEVER...
while True:
    #t0 = time.monotonic_ns()

    #for j in range(100):
    # Get the accelerometer data
    for i in range(fft_size):
        """read x, y, and z as a tuple"""
        acceleration_xyz = accelerometer.acceleration

        samples_x[i] = acceleration_xyz[0]
        samples_y[i] = acceleration_xyz[1]
        samples_z[i] = acceleration_xyz[2]
    #t1 = time.monotonic_ns()
    #r = (t1 - t0) * 1e-6 / 100
    #print("Gathering data takes %8.3fms" % r)
    
    # Suppress DC by subtracting the mean value from the sampled data and apply the Blackman Window
    samples_x = (samples_x - np.mean(samples_x)) * bw
    samples_y = (samples_y - np.mean(samples_y)) * bw
    samples_z = (samples_z - np.mean(samples_z)) * bw

    #t0 = time.monotonic_ns()

    #for i in range(100):
    # Compute the spectral density (voltage mode, i.e. 20*log10(value)) of the sampled data
    """
    spectrogram_x = spectrogram(samples_x)
    spectrogram_x = spectrogram_x[1 : (fft_size // 2) - 1]
    spectrogram_x_dB = 20 * np.log10(spectrogram_x)

    spectrogram_y = spectrogram(samples_y)
    spectrogram_y = spectrogram_y[1 : (fft_size // 2) - 1]
    spectrogram_y_dB = 20 * np.log10(spectrogram_y)
    """
    
    spectrogram_z = spectrogram(samples_z)
    spectrogram_z = spectrogram_z[1 : (fft_size // 2) - 1]
    spectrogram_z_dB = 20 * np.log10(spectrogram_z)

    #t1 = time.monotonic_ns()
    #r = (t1 - t0) * 1.e-6 / 100
    #print("Processing data takes %8.3fms" % r)

    # Print the spectrograms so they can be displayed in the Mu plotter
    for i in range(0, (fft_size // 2) - 2):
        """
        print(
            "(%f, %f, %f)"
            % (spectrogram_x_dB[i], spectrogram_y_dB[i], spectrogram_z_dB[i])
        )
        time.sleep(0.04)
        """
        print("%i," % int(spectrogram_z_dB[i]), end = " ")
    print()


Locked
Please be positive and constructive with your questions and comments.

Return to “Feather - Adafruit's lightweight platform”