10 years, 2 months ago.

mbed-dsp FFT : split magnitude bins?

When I use the code mentioned below (on a FRDM-KL25Z board), I noticed something odd in the magnitudes returned by arm_cmplx_mag_f32:
When i apply a sinusoidal sweep (0..20kHz) to the audio input, i see the max magnitude value go from bin 1 to bin 8 when f = 0 to 10kHz and go backwards from bin 8 to bin 1 when f = 10 to 20kHz.
Unless I mis-interpreted something, the number of magnitude bins is equal to the FFT size?

Assumption, based on my tests
16 magnitude bins are returned when we define a 16 bin FFT, allbeit they are separated into 2 16 bin outputs where only the lower 8 bins are significant? In other words, as the magnitudes are returned in a 16 bin array where the lower half of the array is mirrored into the upper half, the magnitudes output is split into 2 parts:
If frequency @ max magnitude <= samplerate/2, return the lower half of the magnitudes.
If frequency @ max magnitude >= samplerate/2, return the upper half of the magnitudes.
If this is true, how can we determine the corresponding frequency of the returned magnitudes??

Below is a copy of the output taken at 10kHz - Notice the max value for 10kHz is at the 8th bin. If only 8 bins are used, i would expect this to be the 4th bin or, when 16 bins are used, have the upper half magnitudes reflect the 10 to 20kHz range instead of the current mirrored values?

-13.53     0.00     1.09    13.97     4.73    23.37     3.00    16.20   -14.50     0.41   -44.25   -10.98   -94.92   -12.72  -266.36    -3.09  1837.12     0.00  -266.36     3.09   -94.92    12.72   -44.25    10.98   -14.50    -0.41     3.00   -16.20     4.73   -23.37     1.09   -13.97
13.53    14.01    23.84    16.48    14.51    45.59    95.77   266.38  1837.12   266.38    95.77    45.59    14.51    16.48    23.84    14.01
MAX value at magnitudes[8] : +1837.12

#include "mbed.h"
#include <ctype.h>
#include "arm_math.h"
#include "arm_const_structs.h"

Serial pc(USBTX, USBRX);

AnalogIn   left(PTC2);
AnalogIn   right(PTB3);

int SAMPLE_RATE_HZ = 40000;             // Sample rate of the audio in hertz.
const int FFT_SIZE = 16;                // Size of the FFT.

const static arm_cfft_instance_f32 *S;
Ticker samplingTimer;
float samples[FFT_SIZE*2];
float magnitudes[FFT_SIZE];
int sampleCounter = 0;

void samplingCallback()
{
    // Read from the ADC and store the sample data
    samples[sampleCounter] = (1023 * left) - 511.0;
    // Complex FFT functions require a coefficient for the imaginary part of the input.
    // Since we only have real data, set this coefficient to zero.
    samples[sampleCounter+1] = 0.0;
    // Update sample buffer position and stop after the buffer is filled
    sampleCounter += 2;
    if (sampleCounter >= FFT_SIZE*2) {
        samplingTimer.detach();
    }
}

void samplingBegin()
{
    // Reset sample buffer position and start callback at necessary rate.
    sampleCounter = 0;
    samplingTimer.attach_us(&samplingCallback, 1000000/SAMPLE_RATE_HZ);
}

bool samplingIsDone()
{
    return sampleCounter >= FFT_SIZE*2;
}

int main()
{
    // Set up serial port.
    pc.baud (38400);

    // Init arm_ccft_32
    switch (FFT_SIZE)
    {
    case 16:
        S = & arm_cfft_sR_f32_len16;
        break;
    case 32:
        S = & arm_cfft_sR_f32_len32;
        break;
    case 64:
        S = & arm_cfft_sR_f32_len64;
        break;
    case 128:
        S = & arm_cfft_sR_f32_len128;
        break;
    case 256:
        S = & arm_cfft_sR_f32_len256;
        break;
    case 512:
        S = & arm_cfft_sR_f32_len512;
        break;
    case 1024:
        S = & arm_cfft_sR_f32_len1024;
        break;
    case 2048:
        S = & arm_cfft_sR_f32_len2048;
        break;
    case 4096:
        S = & arm_cfft_sR_f32_len4096;
        break;
    }
    float maxValue = 0.0f;
    unsigned int testIndex = 0;

    // Begin sampling audio
    samplingBegin();

    while(1)
    {
        // Calculate FFT if a full sample is available.
        if (samplingIsDone())
        {
            // Run FFT on sample data.
            arm_cfft_f32(S, samples, 0, 1);
            for(int i = 0;i < FFT_SIZE*2;++i)
                printf("% 8.2f ",samples[i]);
            printf("\r\n");
            // Calculate magnitude of complex numbers output by the FFT.
            arm_cmplx_mag_f32(samples, magnitudes, FFT_SIZE);
            for(int i = 0;i < FFT_SIZE;++i)
                printf("% 8.2f ",magnitudes[i]);
            printf("\r\n");
            arm_max_f32(magnitudes, FFT_SIZE, &maxValue, &testIndex);
            printf("MAX value at magnitudes[%d] : %+ 8.2f\r\n", testIndex, maxValue);
            // Wait for user confirmation to restart audio sampling.
            pc.getc();
            samplingBegin();
        }
    }
}

It would be nice to see the raw ADC samples rather than just FFT outputs.

For purely real inputs, you will always have ambiguity in the FFT outputs. That is, you'll always see the FFT magnitude response mirrored around 1/2 the sample frequency. From your post, it looks like you have purely real inputs. It looks to me like your example input signal frequency is not 10 kHz (assuming you're true sample rate is 40 kHz). A 10 kHz signal should show peaks in bins 4 (4 * 40 kHz / 16 = 10 kHz) and 12 (12 * 40 kHz / 16 = 30 kHz) for a 16 point FFT. Note that I always reference bins starting from 0. Without the raw ADC samples it's hard to say what going on with your processing, though.

posted by David G 03 Mar 2014

2 Answers

10 years, 2 months ago.

It depends on the definition, but in principle the top half of a FFT is always the mirror of the bottom half (sometimes it might not be returned since it is a bit pointless). And since an FFT has the same length as its input, with FFT-size of 16 you will have 8 useful bins.

That above 10kHz it goes backwards is Nyquist: http://en.wikipedia.org/wiki/Nyquist%E2%80%93Shannon_sampling_theorem. You need to sample at least two times faster than the highest frequency you have at your input, otherwise it starts folding, which is what you see. Now you are sampling at 40kHz, so that should only happen above 20kHz. That I don't understand yet :). Does the loop indeed runs at the expected speed?

Accepted Answer

Accepted answer - Detailed solution : see my reply below.

posted by Frank Vannieuwkerke 03 Mar 2014
10 years, 2 months ago.

@Erik
Your input solved the 'mystery' The sampling time is way off.
The entire loop takes 804us (measured using Timer), so there's 400us overhead. Exactly twice the time needed to sample (16 x 25us).
The loop time remains at the same value even when the sampling frequency is raised to 80kHz. Theoretically, it should lower by 200us.
I measured one analog read cycle (using Timer) : it takes 40us. I neglected to check the time it takes to get a sample.
So i need to port your FastAnalogIn code to the KL25Z in order to get it working (i don't know if this speed is attainable - need to check the datasheet first).

@David
Yes, i'm only using real input, below is a full acquisition list.
The input signal is 10kHz, generated with this tool : https://code.google.com/p/audiotools/
Eriks remark about the sampling time pointed me into the right direction. I was searching in the wrong places as the real sampling time was, coincidentally, exactly twice the wanted sample time.
ouch, an annoying one - i should have noticed this when i doubled the sampling frequency.

Information

Conclusion : theoretically, the max sampling frequency is 25kHz, but to be on the safe side, we should limit it to 20kHz (code overhead).
Confirmed : a real test @ 20kHz sampling rate works without problems.

loop time :  804us
Raw input  :  -64.34     0.00    69.41     0.00   -75.33     0.00    84.46     0.00   -87.11     0.00    86.02     0.00   -87.35     0.00    84.29     0.00   -82.76     0.00    78.60     0.00   -79.40     0.00    70.52     0.00   -65.26     0.00    57.61     0.00   -52.04     0.00    43.34     0.00
FFT output :  -19.32     0.00    -5.59    -7.34   -10.96    -7.10   -15.25   -12.89    -5.35    -9.04    -3.61   -13.36    21.51   -37.78    98.13   -95.22 -1167.81     0.00    98.13    95.22    21.51    37.78    -3.61    13.36    -5.35     9.04   -15.25    12.89   -10.96     7.10    -5.59     7.34
magnitudes :   19.32     9.22    13.06    19.97    10.50    13.84    43.48   136.73  1167.81   136.73    43.48    13.84    10.50    19.97    13.06     9.22
MAX value at magnitudes[8] : +1167.81


Thanks guys for the fast replies and for pointing me into the right direction.