Explain FFT example (DSP-69)

[tatulea]

Hi, I am looking to do FFT on a signal coming from the accelerometer. I ahve 512 sample readings and I know the frequency. I used FFTW library in python to do the FFT for the signal and it's working. How can I use this library to get the same thing? I had a look over your example, but I don't understand exactly what x1 and x2 are and how I should modify it.

Regards, Gabriel

[dmitry1945]

Hi Gabriel, The esp-dsp library working with complex numbers, that's why the simplest way for you will be just ignore the x2 input and y2 output, and use only x1 and y1. Your code will look so:

    esp_err_t ret;
    ESP_LOGI(TAG, "Start Example.");
    ret = dsps_fft2r_init_fc32(NULL, CONFIG_DSP_MAX_FFT_SIZE);
    if (ret  != ESP_OK)
    {
        ESP_LOGE(TAG, "Not possible to initialize FFT. Error = %i", ret);
        return;
    }
   // Generate hann window
    dsps_wind_hann_f32(wind, N);
    // Generate input signal for x1 A=1 , F=0.1
    dsps_tone_gen_f32(x1, N, 1.0, 0.16,  0);
    // x1 - is your signal that you have to fill
    // Convert two input vectors to one complex vector
    for (int i=0 ; i< N ; i++)
    {
        y_cf[i*2 + 0] = x1[i] * wind[i]; // Real part is your signal multiply with window
        y_cf[i*2 + 1] = 0; // Imag part is 0
    }
    dsps_fft2r_fc32(y_cf, N);
    unsigned int end_b = xthal_get_ccount();
    // Bit reverse 
    dsps_bit_rev_fc32(y_cf, N);
    // Convert one complex vector to two complex vectors
    dsps_cplx2reC_fc32(y_cf, N);

    // y1_cf - is your result in log scale
    // y2_cf - magnitude of your signal in linear scale

    for (int i = 0 ; i < N/2 ; i++) {
        y1_cf[i] = 10 * log10f((y1_cf[i * 2 + 0] * y1_cf[i * 2 + 0] + y1_cf[i * 2 + 1] * y1_cf[i * 2 + 1])/N);
        y2_cf[i] = ((y1_cf[i * 2 + 0] * y1_cf[i * 2 + 0] + y1_cf[i * 2 + 1] * y1_cf[i * 2 + 1])/N);
    }

    // Show power spectrum in 64x10 window from -60 to 0 dB from 0..N/2 samples
    ESP_LOGW(TAG, "Signal x1 in log scale");
    dsps_view(y1_cf, N/2, 64, 10,  -60, 40, '|');
    ESP_LOGW(TAG, "Signal x1 in absolute scale");
    dsps_view(y2_cf, N/2, 64, 10,  0, 2, '|');

[tatulea]

Hi Dmitry,

I am using this xlsx file as input: 512.xlsx It has 512 values and the FFT result that I get using the python script is this: fft_output.xlsx

How can I use you example to get the same result?

[dmitry1945]

Can you show me your python script?

[tatulea]

The code that I am using with python is from mide.com

It is not too clean and I've modified few things to accommodate my data input.

import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import tkinter as tk
from tkinter import filedialog
import time
import pyfftw
import os

#####
# Functions for testing various FFT methods
#####

def fftw_builder_test(input_data):
    pyfftw.forget_wisdom()  # This is just here to keep the tests honest, normally pyfftw will remember setup parameters and go more quickly when it is run a second time
    a = np.array(input_data,
                 dtype='float32')  # This turns the input list into a numpy array. We can make it a 32 bit float because the data is all real, no imaginary component
    fft_obj = pyfftw.builders.rfft(a)  # This creates an object which generates the FFT
    return fft_obj()  # And calling the object returns the FFT

def fftw_fast_builder_test(input_data):  # See fftw_builder_test for comments
    pyfftw.forget_wisdom()
    a = np.array(input_data, dtype='float32')
    fft_obj = pyfftw.builders.rfft(a,
                                   planner_effort='FFTW_ESTIMATE')  # FFTW_ESTIMATE is a lower effort planner than the default. This seems to work more quickly for over 8000 points
    return fft_obj()

def fftw_test(input_data):
    pyfftw.forget_wisdom()  # This is just here to keep the tests honest
    outLength = len(
        input_data) // 2 + 1  # For a real FFT, the output is symetrical. fftw returns only half the data in this case
    a = pyfftw.empty_aligned(len(input_data),
                             dtype='float32')  # This is the input array. It will be cleared when the fft object is created
    outData = pyfftw.empty_aligned(outLength,
                                   dtype='complex64')  # This is the output array. Not that the size and type must be appropriate
    fft_obj = pyfftw.FFTW(a, outData, flags=('FFTW_ESTIMATE',),
                          planning_timelimit=1.0)  # NB: The flags tuple has a comma
    a[:] = np.array(input_data,
                    dtype='float32')  # We have to fill the array after fft_obj is created. NB: 'a[:] =' puts data into the existing array, 'a =' creates a new array
    return fft_obj()  # Calling the object returns the FFT of the data now in a. The result is also in outData

def fftw_test_complex(
        input_data):  # This is fftw_test but running a complex FFT as opposed to a real input FFT. See fftw_test for comments
    pyfftw.forget_wisdom()
    outLengthMod = len(input_data) // 2 + 1  # Size of expected return data
    outLength = len(input_data)  # For a complex input FFT, we get more data
    a = pyfftw.empty_aligned(len(input_data),
                             dtype='complex64')  # The FFTW determines the type of FFT by the type of input
    outData = pyfftw.empty_aligned(outLength, dtype='complex64')
    fft_obj = pyfftw.FFTW(a, outData, flags=('FFTW_ESTIMATE',), planning_timelimit=1.0)
    a[:] = np.array(input_data, dtype='complex64')
    return fft_obj()[0:outLengthMod]

def fftw_test_default(input_data):  # This is fftw_test with different FFTW options. See fftw_test for comments
    pyfftw.forget_wisdom()
    outLength = len(input_data) // 2 + 1
    a = pyfftw.empty_aligned(len(input_data), dtype='float32')
    outData = pyfftw.empty_aligned(outLength, dtype='complex64')
    fft_obj = pyfftw.FFTW(a, outData)
    a[:] = np.array(input_data, dtype='float32')
    return fft_obj()

def fftw_test_no_limit(input_data):  # This is fftw_test with different FFTW options. See fftw_test for comments
    pyfftw.forget_wisdom()
    outLength = len(input_data) // 2 + 1
    a = pyfftw.empty_aligned(len(input_data), dtype='float32')
    outData = pyfftw.empty_aligned(outLength, dtype='complex64')
    fft_obj = pyfftw.FFTW(a, outData, flags=('FFTW_ESTIMATE',))
    a[:] = np.array(input_data, dtype='float32')
    return fft_obj()

def fftw_test_no_flag(input_data):  # This is fftw_test with different FFTW options. See fftw_test for comments
    pyfftw.forget_wisdom()
    outLength = len(input_data) // 2 + 1
    a = pyfftw.empty_aligned(len(input_data), dtype='float32')
    outData = pyfftw.empty_aligned(outLength, dtype='complex64')
    fft_obj = pyfftw.FFTW(a, outData, planning_timelimit=1.0)
    a[:] = np.array(input_data, dtype='float32')
    return fft_obj()

def numpy_test(inputData):
    return np.fft.rfft(inputData)  # Numpy is nice and simple!

#####
# Helper Functions
#####

def write_csv(csv_file_name, headers, data):  # Helper function for writing to the CSV file
    outString = ""
    for header in headers:
        outString += "%s," % (data.get(header, ""))
    outString += "\n"
    with open(csv_file_name, 'a') as f:
        f.write(outString)

# Prompt user for file
###
# RUN TIME OPTION!
# Set include_plots to True if you want to actually see the plot outputs. They will all be stored up and displayed
# at the end of the test
include_plots = True
if include_plots:
    root = tk.Tk()
    root.withdraw()
file_paths = filedialog.askopenfilenames(filetypes=[("Two Column CSV", "*.csv")])
out_path = filedialog.asksaveasfilename(defaultextension=".csv")
print(file_paths)

# CSV Setup
output_headers = ["File", "File Size", "Number of Points", "FFT Time", "FFT Function", "Load Time", "Plot Time",
                  "RMS Time"]
with open(out_path, 'w+') as f:
    f.write(",".join(output_headers) + "\n")

# Tests to run:
# Format { "Name" : Function, etc } Test order is undefined
# test_functions = {"FFTW" : fftw_test, "FFTW-Complex" : fftw_test_complex, "FFTW-Builder": fftw_builder_test, "FFTW-Default" : fftw_test_default, "FFTW-NoLimit": fftw_test_no_limit, "FFTW-NoFlag": fftw_test_no_flag}
test_functions = {"FFTW": fftw_test, "Builder": fftw_fast_builder_test}
figure_count = 1
# Loop through each file
for file_path in file_paths:
    # output_data holds the test info and results. This is written to the file after the FFT is complete
    output_data = {"File": os.path.split(file_path)[-1]}
    output_data["File Size"] = os.path.getsize(file_path)
    tic = time.clock()
    df = pd.read_csv(file_path, delimiter='\t', header=None, names=["x"])
    x = df["x"]
    toc = time.clock()
    output_data["Load Time"] = toc - tic
    print("Load Time:", toc - tic)

    N = np.int(np.prod(x.shape))  # length of the array
    # Determine variables
    Fs = 3333
    T = 1 / Fs;
    output_data["Number of Points"] = N
    print("# Samples:", N)

    # Compute and Plot FFT
    # Run through each FFT function in test_functions, and test each one.
    # File data is not reloaded
    fftData = None
    xf = None
    for func_name in test_functions:
        output_data["FFT Function"] = func_name
        tic = time.clock()
        fftData = test_functions[func_name](x)  # Execute the FFT
        xf = np.linspace(0.0, 1.0 / (2.0 * T), len(fftData))
        if include_plots:
            plt.figure(figure_count)
            figure_count += 1
            plt.plot(xf, 2.0 / N * np.abs(fftData))  # [0:np.int(N/2)]
            plt.grid()
            plt.xlabel('Frequency (Hz)')
            plt.ylabel('Accel (g)')
            plt.title('FFT (%s) - %s' % (func_name, file_path))
        toc = time.clock()
        output_data["FFT Time"] = toc - tic
        print("FFT Time:", toc - tic)
        write_csv(out_path, output_headers, output_data)

    output_data = []
    for i, elem in enumerate(xf):
        output_data.append(str(elem) + "," + str(2.0 / N * np.abs(fftData[i])))
    with open('/home/gabi/iesite.csv', 'w') as f:
        f.write("\n".join(output_data))
if include_plots:
    print("Plotting...")
    plt.show()
print("Done!")

On lines 176, 177 I get the output to csv and I think that is what you are looking for

[dmitry1945]

OK. Just take y2_cf output and not use window from example that I have included. Like this:

    esp_err_t ret;
    ESP_LOGI(TAG, "Start Example.");
    ret = dsps_fft2r_init_fc32(NULL, CONFIG_DSP_MAX_FFT_SIZE);
    if (ret  != ESP_OK)
    {
        ESP_LOGE(TAG, "Not possible to initialize FFT. Error = %i", ret);
        return;
    }
   // Generate hann window
    // Generate input signal for x1 A=1 , F=0.1
    dsps_tone_gen_f32(x1, N, 1.0, 0.16,  0);
    // x1 - is your signal that you have to fill
    // Convert two input vectors to one complex vector
    for (int i=0 ; i< N ; i++)
    {
        y_cf[i*2 + 0] = x1[i]; // Real part is your signal
        y_cf[i*2 + 1] = 0; // Imag part is 0
    }
    dsps_fft2r_fc32(y_cf, N);
    unsigned int end_b = xthal_get_ccount();
    // Bit reverse 
    dsps_bit_rev_fc32(y_cf, N);
    // Convert one complex vector to two complex vectors
    dsps_cplx2reC_fc32(y_cf, N);

    // y2_cf - magnitude of your signal in linear scale

    for (int i = 0 ; i < N/2 ; i++) {
       y2_cf[i] = ((y_cf[i * 2 + 0] * y_cf[i * 2 + 0] + y_cf[i * 2 + 1] * y_cf[i * 2 + 1])/N);
    }

    // Show power spectrum in 64x10 window from -60 to 0 dB from 0..N/2 samples
    ESP_LOGW(TAG, "Signal x1 in absolute scale");
    dsps_view(y2_cf, N/2, 64, 10,  0, 2, '|');

[tatulea]

I am printing y2_cf elements, but they don't look like the ones from xlsx file. Shouldn't they be the same?

[dmitry1945]

They should. Can you show them here?

[dmitry1945]

And, the last example is different then the firs.

[tatulea]

The code that I am using is:

    static float x[512] = {1.6,1.3,1,0.6,0.3,0.1,0,0.1,0.2,0.5,0.8,1.2,1.5,1.8,2,2,2,1.8,1.5,1.1,0.8,0.5,0.2,0.1,0,0.1,0.4,0.7,1,1.4,1.7,1.9,2,2,1.9,1.6,1.3,0.9,0.6,0.3,0.1,0,0.1,0.3,0.5,0.9,1.2,1.5,1.8,2,2,1.9,1.7,1.4,1.1,0.7,0.4,0.2,0.1,0,0.2,0.4,0.7,1.1,1.4,1.7,1.9,2,2,1.8,1.6,1.2,0.9,0.5,0.3,0.1,0,0.1,0.3,0.6,0.9,1.3,1.6,1.8,2,2,1.9,1.7,1.4,1,0.7,0.4,0.2,0,0.1,0.2,0.4,0.8,1.1,1.5,1.7,1.9,2,2,1.8,1.5,1.2,0.8,0.5,0.2,0.1,0,0.1,0.3,0.6,1,1.3,1.6,1.9,2,2,1.9,1.6,1.3,1,0.6,0.3,0.1,0,0.1,0.2,0.5,0.8,1.2,1.5,1.8,2,2,1.9,1.7,1.5,1.1,0.8,0.5,0.2,0.1,0,0.1,0.4,0.7,1,1.4,1.7,1.9,2,2,1.8,1.6,1.3,0.9,0.6,0.3,0.1,0,0.1,0.3,0.5,0.9,1.2,1.6,1.8,2,2,1.9,1.7,1.4,1.1,0.7,0.4,0.2,0,0,0.2,0.4,0.7,1.1,1.4,1.7,1.9,2,2,1.8,1.6,1.2,0.9,0.5,0.3,0.1,0,0.1,0.3,0.6,0.9,1.3,1.6,1.8,2,2,1.9,1.7,1.4,1,0.7,0.4,0.1,0,0.1,0.2,0.5,0.8,1.1,1.5,1.8,1.9,2,2,1.8,1.5,1.2,0.8,0.5,0.2,0.1,0,0.1,0.3,0.6,1,1.3,1.6,1.9,2,2,1.9,1.6,1.3,1,0.6,0.3,0.1,0,0.1,0.2,0.5,0.8,1.2,1.5,1.8,2,2,1.9,1.7,1.5,1.1,0.8,0.4,0.2,0.1,0,0.1,0.4,0.7,1,1.4,1.7,1.9,2,2,1.8,1.6,1.3,0.9,0.6,0.3,0.1,0,0.1,0.3,0.5,0.9,1.2,1.6,1.8,2,2,1.9,1.7,1.4,1.1,0.7,0.4,0.2,0,0,0.2,0.4,0.7,1.1,1.4,1.7,1.9,2,2,1.8,1.5,1.2,0.9,0.5,0.3,0.1,0,0.1,0.3,0.6,0.9,1.3,1.6,1.9,2,2,1.9,1.7,1.3,1,0.7,0.4,0.1,0,0.1,0.2,0.5,0.8,1.1,1.5,1.8,2,2,2,1.8,1.5,1.2,0.8,0.5,0.2,0.1,0,0.1,0.3,0.6,1,1.4,1.7,1.9,2,2,1.9,1.6,1.3,0.9,0.6,0.3,0.1,0,0.1,0.2,0.5,0.8,1.2,1.5,1.8,2,2,1.9,1.7,1.4,1.1,0.7,0.4,0.2,0,0,0.2,0.4,0.7,1.1,1.4,1.7,1.9,2,2,1.8,1.6,1.2,0.9,0.6,0.3,0.1,0,0.1,0.3,0.6,0.9,1.3,1.6,1.8,2,2,1.9,1.7,1.4,1,0.7,0.4,0.2,0,0.1,0.2,0.4,0.8,1.1,1.5,1.7,1.9,2,2,1.8,1.5,1.2,0.8,0.5,0.2,0.1,0,0.1,0.3,0.6,1,1.3,1.6,1.9,2,2,1.9,1.6,1.3,1,0.6,0.3,0.1,0,1.7,1.9,2,2,1.8,1.6,1.3,0.9,0.6,0.3,0.1,0,0.1,0.3,0.5,0.9,1.2,1.6,1.8,2,2,1.9,1.7,1.4,1.1,0.7,0.4,0.2,0,0,0.2,0.4,0.7,1.1,1.4,1.7,1.9,2,2,1.8,1.5,1.2,0.9,0.5,0.3,0.1,0,0.1,0.3,0.6,0.9,1.3};

    esp_err_t ret;
    ESP_LOGI(TAG, "Start Example.");
    ret = dsps_fft2r_init_fc32(NULL, CONFIG_DSP_MAX_FFT_SIZE);
    if (ret  != ESP_OK)
    {
        ESP_LOGE(TAG, "Not possible to initialize FFT. Error = %i", ret);
        return;
    }
   // Generate hann window
    // Generate input signal for x1 A=1 , F=0.1
    dsps_tone_gen_f32(x, N, 1.0, 0.16,  0);
    // x1 - is your signal that you have to fill
    // Convert two input vectors to one complex vector
    for (int i=0 ; i< N ; i++)
    {
        y_cf[i*2 + 0] = x[i]; // Real part is your signal
        y_cf[i*2 + 1] = 0; // Imag part is 0
    }
    dsps_fft2r_fc32(y_cf, N);
    unsigned int end_b = xthal_get_ccount();
    // Bit reverse 
    dsps_bit_rev_fc32(y_cf, N);
    // Convert one complex vector to two complex vectors
    dsps_cplx2reC_fc32(y_cf, N);

    // y2_cf - magnitude of your signal in linear scale

    for (int i = 0 ; i < N/2 ; i++) {
       y2_cf[i] = ((y_cf[i * 2 + 0] * y_cf[i * 2 + 0] + y_cf[i * 2 + 1] * y_cf[i * 2 + 1])/N);
       printf("%f\r\n", y2_cf[i]);
    }

    // Show power spectrum in 64x10 window from -60 to 0 dB from 0..N/2 samples
    ESP_LOGW(TAG, "Signal x1 in absolute scale");
    dsps_view(y2_cf, N/2, 64, 10,  0, 2, '|');

The output that I am getting is:

0.000244
0.000976
0.000977
0.000979
0.000982
0.000986
0.000991
0.000996
0.001003
0.001010
0.001019
0.001028
0.001038
0.001049
0.001062
0.001075
0.001090
0.001106
0.001123
0.001141
0.001161
0.001182
0.001205
0.001229
0.001256
0.001284
0.001314
0.001346
0.001380
0.001416
0.001456
0.001497
0.001542
0.001590
0.001642
0.001697
0.001756
0.001819
0.001888
0.001961
0.002040
0.002126
0.002218
0.002318
0.002426
0.002543
0.002671
0.002811
0.002963
0.003130
0.003313
0.003514
0.003737
0.003983
0.004257
0.004562
0.004904
0.005288
0.005722
0.006215
0.006777
0.007422
0.008168
0.009036
0.010053
0.011258
0.012697
0.014437
0.016567
0.019212
0.022552
0.026853
0.032523
0.040209
0.050992
0.066789
0.091265
0.132163
0.208265
0.375517
0.869075
3.787957
501.346039
2.753774
0.743241
0.339379
0.193663
0.125103
0.087471
0.064615
0.049699
0.039429
0.032057
0.026587
0.022416
0.019163
0.016577
0.014487
0.012773
0.011352
0.010158
0.009147
0.008282
0.007538
0.006892
0.006327
0.005831
0.005393
0.005005
0.004658
0.004347
0.004068
0.003816
0.003588
0.003380
0.003191
0.003019
0.002861
0.002715
0.002581
0.002458
0.002344
0.002238
0.002140
0.002048
0.001963
0.001883
0.001809
0.001740
0.001674
0.001613
0.001556
0.001501
0.001450
0.001402
0.001357
0.001314
0.001273
0.001234
0.001198
0.001163
0.001130
0.001099
0.001069
0.001041
0.001013
0.000988
0.000963
0.000939
0.000917
0.000895
0.000875
0.000855
0.000836
0.000818
0.000800
0.000784
0.000768
0.000752
0.000737
0.000723
0.000710
0.000696
0.000684
0.000671
0.000660
0.000648
0.000637
0.000627
0.000617
0.000607
0.000597
0.000588
0.000579
0.000571
0.000562
0.000554
0.000546
0.000539
0.000532
0.000525
0.000518
0.000511
0.000505
0.000499
0.000493
0.000487
0.000481
0.000476
0.000471
0.000465
0.000461
0.000456
0.000451
0.000447
0.000442
0.000438
0.000434
0.000430
0.000426
0.000422
0.000418
0.000415
0.000411
0.000408
0.000405
0.000402
0.000399
0.000396
0.000393
0.000390
0.000388
0.000385
0.000382
0.000380
0.000378
0.000375
0.000373
0.000371
0.000369
0.000367
0.000365
0.000363
0.000362
0.000360
0.000358
0.000357
0.000355
0.000354
0.000352
0.000351
0.000350
0.000348
0.000347
0.000346
0.000345
0.000344
0.000343
0.000342
0.000341
0.000340
0.000340
0.000339
0.000338
0.000338
0.000337
0.000337
0.000336
0.000336
0.000336
0.000335
0.000335
0.000335
0.000335
0.000334
0.000334
W (558) main: Signal x1 in absolute scale
I (568) view: Data min[0] = 0.000244, Data max[82] = 501.346039
 ________________________________________________________________
0                    |                                           |
1                    |                                           |
2                    |                                           |
3                    |                                           |
4                    |                                           |
5                    ||                                          |
6                     |                                          |
7                   | |                                          |
8||||||||||||||||||||  |||||||||||||||||||||||||||||||||||||||||||
9                                                                |
 0123456789012345678901234567890123456789012345678901234567890123

[dmitry1945]

Yes! :) y2_cf[i] = sqrt((y_cf[i 2 + 0] y_cf[i 2 + 0] + y_cf[i 2 + 1] y_cf[i 2 + 1]))/N;

))))

[tatulea]

Is this working for you? The first values from my expected output are

x | y
-- | --
6.50976563 | 0.01152959
13.0195313 | 0.01720576
19.5292969 | 0.0210918
26.0390625 | 0.02718033
32.5488281 | 0.0254768
39.0585938 | 0.02137978
45.5683594 | 0.01765214
52.078125 | 0.01952043
58.5878906 | 0.00869537
65.0976563 | 0.01339019
71.6074219 | 0.0224928
78.1171875 | 0.02958141

And the values that I get from ESP are

0.000690
0.001381
0.001382
0.001383
0.001385
0.001388
0.001391
0.001395
0.001400
0.001405
0.001410
0.001417

[tatulea]

I am not too much into math and this is why I don't know exactly how to solve this. My final point should be the fact that the x input that I am using should return me a peak frequency of about 200MHz

[dmitry1945]

dsps_tone_gen_f32(x, N, 1.0, 0.16, 0); is an error. You overwrite your input signal with tone. Remove this line.

[tatulea]

I tried, but I still got the wrong values

[this-username-is-taken]

Hi Gabriel, The esp-dsp library working with complex numbers, that's why the simplest way for you will be just ignore the x2 input and y2 output, and use only x1 and y1. Your code will look so:

    esp_err_t ret;
    ESP_LOGI(TAG, "Start Example.");
    ret = dsps_fft2r_init_fc32(NULL, CONFIG_DSP_MAX_FFT_SIZE);
    if (ret  != ESP_OK)
    {
        ESP_LOGE(TAG, "Not possible to initialize FFT. Error = %i", ret);
        return;
    }
   // Generate hann window
    dsps_wind_hann_f32(wind, N);
    // Generate input signal for x1 A=1 , F=0.1
    dsps_tone_gen_f32(x1, N, 1.0, 0.16,  0);
    // x1 - is your signal that you have to fill
    // Convert two input vectors to one complex vector
    for (int i=0 ; i< N ; i++)
    {
        y_cf[i*2 + 0] = x1[i] * wind[i]; // Real part is your signal multiply with window
        y_cf[i*2 + 1] = 0; // Imag part is 0
    }
    dsps_fft2r_fc32(y_cf, N);
    unsigned int end_b = xthal_get_ccount();
    // Bit reverse 
    dsps_bit_rev_fc32(y_cf, N);
    // Convert one complex vector to two complex vectors
    dsps_cplx2reC_fc32(y_cf, N);

    // y1_cf - is your result in log scale
    // y2_cf - magnitude of your signal in linear scale

    for (int i = 0 ; i < N/2 ; i++) {
        y1_cf[i] = 10 * log10f((y1_cf[i * 2 + 0] * y1_cf[i * 2 + 0] + y1_cf[i * 2 + 1] * y1_cf[i * 2 + 1])/N);
        y2_cf[i] = ((y1_cf[i * 2 + 0] * y1_cf[i * 2 + 0] + y1_cf[i * 2 + 1] * y1_cf[i * 2 + 1])/N);
    }

    // Show power spectrum in 64x10 window from -60 to 0 dB from 0..N/2 samples
    ESP_LOGW(TAG, "Signal x1 in log scale");
    dsps_view(y1_cf, N/2, 64, 10,  -60, 40, '|');
    ESP_LOGW(TAG, "Signal x1 in absolute scale");
    dsps_view(y2_cf, N/2, 64, 10,  0, 2, '|');

This is great! This should be in the documentary. Made things so much easier to understand.

Btw, why is there not a real-input-only implementation of FFT?

[nakanaela]

From my experience with rFFT algorithms in the CMSIS DSP library, there really is no such thing as a 'real' only fft. Those functions do a bit of hand-waving to allow you to use a vector of real only inputs and creates a secondary vector of 2sample size with every i2+1 index set to zero to account for the complex portion. The benefit is that you don't have to do it yourself, but the resulting firmware files were easily twice what I get from zeroing the complex index myself and directly using the corresponding cFFT function.

If you are sampling real data, it is pretty easy to account for the complex value being '0' and in my experience has been faster than refactoring a vector of real inputs to a vector of real+cmplx values. If you initialize your input array to '0' from the start, then all you have to worry about is adding your real inputs at every [2sample#] as [2sample#+1] should already be '0'.

[drewandre]

Hi @dmitry1945 - I was hoping you could elaborate on the fft example a bit more in the context of audio streaming. I have the example referred to above working great but I'd like to use the audio stream from an esp-adf i2s stream and I don't see any examples in any esp library nor have I seen any helpful forum posts which is surprising. I imagine others would want to access the i2s buffer to do some audio visual processing or similar.

At a high level, I'm trying to stream audio to the esp32 lyrat 4.3 board, run an fft, and display the spectrum across an led strip. I have each element working in isolation (the bluetooth audio stream, writing i2s data to the lyrat audio codec, outputting audio from the codec, running non-audio reactive led animations, and running the above fft example on a single tone using dsps_tone_gen_f32).

I've tried a few approaches including 1) trying to read the i2s raw byte array using i2s_read and iterating over that, 2) creating an audio pipeline and passing the byte array over to the fft function, and 3) implementing an i2s write callback.

Each of these either doesn't compile or work, but I'm also not sure if there's something obvious I'm missing or if I'm not on the right track at all. Any pointers here on what methods I could leverage from the esp-adf library to pass the i2s buffer over to the dsps_fft2r_fc32 function?

Thanks!

[dmitry1945]

Hi @drewandre ,

If I understood correct, you have next application: I2S_read/bluetooth read, fft, some adf codec, I2S_write. I would suggest to try next: prepare the buffer with dsps_tone_gen_f32 and pass this buffer as an predefined input to the dsps_fft2r_fc32 inside your application instead of buffer with real data. In this case you should see the expected result. If you see the expected result, then we know, that the problem is not in the application and in functions calls, but in the data format. If the problem in data format, then just convert your input to the float and use or use dsps_fft2r_fc32 as in example, or use dsps_fft2r_sc16 instead. If you don't see the expected result with predefined data, then the problem in the application and you have to care about execution flow. Please let me know if it so, I will provide suggestions.

Regards, Dmitry

[drewandre]

Thanks for the quick response @dmitry1945. You are correct that the application is:

                          |-> FFT from i2s stream -> LED animations from FFT data
bluetooth audio stream -> |
                          |-> i2s write -> audio codec headphone output

I have already been using dsps_tone_gen_f32 as you suggested and it works well. I was using this just to make sure the fft logic is running correctly.

This question is actually more of an esp-adf library question so I opened an issue on the esp-adf page. But I'm wondering if you have any pointers on how to pass an i2s stream to dsps_fft2r_fc32. I assumed this library would play a little nicer with esp-adf like the teensy audio library which has an i2s audio element that can be "patched" directly to the fft but it's clear that these are two completely separate libraries. For example, the fft functions in esp-dsp accept floating or fixed buffers, but the audio element in esp-adf works with char buffers.

The repo I'm working on is basically the a2dp_gatts_coex esp-idf example. I won't be using the fft for audio processing - I would be using it to control LED animations so the fft wouldn't have to be part of the audio pipeline.

So my question is do you have any tips on how to access the char buffer array from the esp-adf audio element so I can use it with the esp-dsp fft? When I get this working I would be happy to create a PR with an example on ways to use this library with esp-adf. If you feel this question isn't appropriate for this library, feel free to close it this. I would appreciate any tips!

[dmitry1945]

Hi @drewandre ,

I think the best way to do FFT is to run it on separate thread. Make thread with semaphore, and when data comes to the I2s stream callback just copy the i2s buffer to the fft processing buffer, and then trigger the semaphore if the FFT thread is not busy. Something like this:

bool fft_busy = false;
byte data_for_fft[NNN];
void some_i2s_callback()
{
     ... fill buffer with data for fft - data_for_fft
     if (false == fft_busy) xSemaphoreGive(fft_task_sem);
     .... some other processing
}
void FFT_thread()
    while (1)
    {
        fft_busy  = false;
        xSemaphoreTake(fft_task_sem, 10000000);
        fft_busy = true;
        .... here convert i2s buffer (data_for_fft) to fft complex float
        .... process buffer in fft and calculate magnitude
        .... show buffer
    }
}

Regards, Dmitry

[drewandre]

EDIT: I may have solved this by referencing the logic in this hackaday post. Testing now.

EDIT 2: I've solved my issue after stumbling across that hackaday post. I've changed the loop in step 5 to:

int t = 0;
static int16_t sample_l_int = 0;
static int16_t sample_r_int = 0;
static float sample_l_float = 0.0f;
static float sample_r_float = 0.0f;
static float in = 0.0f;

for (uint32_t i = 0; i < N_SAMPLES; i += 4)
{
  sample_l_int = (int16_t)((*(fft_buffer + i + 1) << 8) | *(fft_buffer + i));
  sample_r_int = (int16_t)((*(fft_buffer + i + 3) << 8) | *(fft_buffer + i + 2));
  sample_l_float = (float)sample_l_int / 0x8000;
  sample_r_float = (float)sample_r_int / 0x8000;
  in = (sample_l_float + sample_r_float) / 2.0f;
  x1[t] = in * wind[t];
  t++;
}

I've also decreased my window array length by a factor of 4.

Hi @dmitry1945, thanks again for your help. I think I'm pretty close to getting esp-dsp's FFT's running against a bluetooth audio stream (interleaved PCM stereo audio) but dsps_view's serial output just doesn't look right even with a pure 440hz sine wave. I'm still running ESP-IDF 3.3.1 and ESP-ADF v2.0-beta3.

The serial output of dsps_view can be viewed here, and the float buffer passed to dsps_fft2r_fc32 and related functions can be viewed here.

For reference, here is dsps_view's serial output when swapping out my PCM buffer with the one generated by dsps_tone_gen_f32(x1, N, 1.0, 0.16, 0) (only copying this here to show that the FFT logic is correct - it just seems to be a data format issue with my PCM buffer). I do notice here that the the min property printed by dsps_view is -91.273849 for my PCM 440hz sine tone, and -inf when using dsps_tone_gen_f32's output. Unsure if this is relevant.

Do you mind checking my fft logic to see if anything jumps out? fft_controller.c in full can be viewed on my github repo here. At a high level, here is what I'm doing:

1. Initialize an a2dp data callback to copy the outgoing buffer over to fft_buffer to be used for FFT analysis:

   
   #define N_SAMPLES 2048
   uint8_t fft_buffer[N_SAMPLES * 2] // 4096 because that's the length of the esp-adf's a2dp callback buffer

// data is a pointer to a uint8_t array of PCM interleaved stereo audio data // len is always 4096 here void copy_a2dp_buffer_to_fft_buffer(const uint8_t *data, uint32_t len) { if (!fft_running) { memcpy(fft_buffer, data, len); } }

2. Initialize fft, dynamically allocate some buffers, and fill `wind` and `x1` with window and tone data (`dsps_tone_gen_f32` doesn't really need to be called, just added it to init a buffer of data that we know is correct). **I assume dynamically allocating these arrays as well as passing `NULL` as a table_buff parameter is ok? I also assume use of calloc over malloc wouldn't be an issue.**

float x1, wind, y_cf, y1_cf;

void init_fft() { esp_err_t ret = dsps_fft2r_init_fc32(NULL, CONFIG_DSP_MAX_FFT_SIZE); x1 = (float )calloc(N_SAMPLES, sizeof(float)); wind = (float )calloc(N_SAMPLES, sizeof(float)); y_cf = (float )calloc(N_SAMPLES 2, sizeof(float)); y1_cf = &y_cf[0];

dsps_wind_hann_f32(wind, N_SAMPLES); dsps_tone_gen_f32(x1, N_SAMPLES, 1.0, 0.16, 0); }

3. Attach the fft task (I haven't been using semaphores yet):

xTaskCreatePinnedToCore(&calculate_fft, "calculate_fft", 4000, NULL, 5, NULL, 1);

4. Establish an a2dp connection and begin streaming audio
5. In `calculate_fft` which runs every 15 seconds, fill `x1` float array with formatted values from `fft_buffer`.  The simple `map` function used can be viewed [here](https://gist.github.com/drewandre/0bfcb76e5ad785308ea5ec740d857bf8). **This is likely where the issue is since this is where the data is actually mapped to the format `dsps_fft2r_fc32` expects:**

fft_running = true; // set flag which stops the memcpy in step 1

for (int i = 0; i < N_SAMPLES; i++) { // i 2 since stereo interleaved buffer is structured [left][right][left][right]... // so in this example we are only building an fft buffer for the left channel. // we then divide the sample by 255.0 to reduce the range to 0 - 1. float sample = (float)fft_buffer[i 2] / 255.0;

// here we map the sample from 0 - 1 to -1 to 1. x1[i] = map(sample, 0, 1, -1, 1); }

5. Run the fft (same exact logic as the [fft example](https://github.com/espressif/esp-dsp/tree/master/examples/fft) but with a larger `dsps_view` output):

// Convert two input vectors to one complex vector for (int i = 0; i < N_SAMPLES; i++) { y_cf[i 2 + 0] = x1[i] wind[i]; y_cf[i * 2 + 1] = 0; }

dsps_fft2r_fc32(y_cf, N_SAMPLES); // Bit reverse dsps_bit_rev_fc32(y_cf, N_SAMPLES); // Convert one complex vector to two complex vectors dsps_cplx2reC_fc32(y_cf, N_SAMPLES);

for (int i = 0; i < N_SAMPLES / 2; i++) { y1_cf[i] = 10 log10f((y1_cf[i 2 + 0] y1_cf[i 2 + 0] + y1_cf[i 2 + 1] y1_cf[i * 2 + 1]) / N_SAMPLES); }

dsps_view(y1_cf, N_SAMPLES / 2, 128, 20, -60, 40, '|');

fft_running = false; // set flag to allow memcpy to copy PCM buffer to fft_buffer in step 1

vTaskDelay(15000 / portTICK_PERIOD_MS); // only run every 15 sec for debugging

Anything look incorrect to you? I understand this may be more of an issue with how I'm handling the PCM buffer and not an issue with the library, so feel free to just say you're not sure and I'll go ahead and post on some forums. I see in [this forum comment](https://esp32.com/viewtopic.php?t=11145#p45369) that the a2dp callback buffer is "Signed 16 bits, big-edian, 2 channel" for what it's worth.

It would be great to figure this out so other users could more easily get esp-adf and esp-dsp to play nicely together :) Thanks in advance!

[dmitry1945]

Hi @drewandre, Just to be sure, how do you initialize the y1_cf value? I think in our case it's possible to use next: y1_cf[i] = 10 * log10f((y_cf[i * 2 + 0] * y_cf[i * 2 + 0] + y_cf[i * 2 + 1] * y_cf[i * 2 + 1]) / N_SAMPLES); UPD: I've found y1_cf . :)

In this case you can try to remove dsps_cplx2reC_fc32(y_cf, N_SAMPLES); from the code.

dsps_fft2r_fc32(y_cf, N_SAMPLES);
// Bit reverse
dsps_bit_rev_fc32(y_cf, N_SAMPLES);

for (int i = 0; i < N_SAMPLES / 2; i++)
{
  y1_cf[i] = 10 * log10f((y1_cf[i * 2 + 0] * y1_cf[i * 2 + 0] + y1_cf[i * 2 + 1] * y1_cf[i * 2 + 1]) / N_SAMPLES);
}

dsps_view(y1_cf, N_SAMPLES / 2, 128, 20, -60, 40, '|');

[drewandre]

@dmitry1945 Thanks for the tip - I've made this adjustment.

I have a couple last questions for you if you don't mind:

1. Regarding memory management, would you suggest passing my x1 buffer to dsps_fft2r_init_fc32 rather than having it internally allocated?
2. I'm ultimately looking to use all 512 frequency bins to drive LED animations with each bin's amplitude ranging from 0 to 1.0. How would you recommend mapping the output of dsps_bit_rev_fc32 to achieve this?

[dmitry1945]

Hi @drewandre , I will try to answer,

1. The semaphore hold the FFT task until event will come. It just synchronize your BT callback and FFT calculation.
2. The dsps_bit_rev_fc32 is not just a loop. It's doing some more complicated work. For example if length is 256, bit reverse of 10000000b is 00000001b, 00000010b is 01000000b, 00000110b is 01100000b, and so on.
3. I think better to have as it is.

Regards, Dmitry

[wired8]

Hi @drewandre did you ever get this working? I'm looking to calculate 1/3 octaves from a streaming source.

[drewandre]

Yes, it has been over a year since I've worked on the project but I've attached code samples that may be helpful.

audio_analysis_controller.cpp: https://gist.github.com/dandre-hound/952127e65cb5565e3928c5121977132d

audio_analysis_controller.hpp: https://gist.github.com/dandre-hound/e4357b422a144fb8f4ba69ee0a8bdf2a

There’s a lot of extra code in here unrelated to what you’re trying to do, but to calculate octave bands, I created an initial octave_bands array, and then used the pre-defined frequency bounds to group the outputs of the fft analysis. Relevant functions in these files include average_into_octave_bands, calcAvg, and freqToIndex.

[wired8]

Thanks for that @drewandre much appreciated.