search

espressif / esp-dsp

DSP library for ESP-IDF
Apache License 2.0
466 stars 87 forks source link

How do I know the frequencies displayed by 'dsps_view' function? (DSP-140) #88

Closed marco711 closed 1 month ago

marco711 commented 1 month ago

Answers checklist.

General issue report

I'm testing the fft example provided by Espressif. Nevertheless, in my application I'm injecting the audio signal captured by a microphone connected to the ADC input. What I'd like to know is how to read the frequencies displayed by the dsps_view function. With the example provided by Espressif I get:

 ________________________________________________________________
0                                                                |
1                    |                                           |
2                    |                                           |
3                    |    |                                      |
4                    |    |                                      |
5                    |    |                                      |
6                   | |   |                                      |
7                   | |   ||                                     |
8                  || || | |                                     |
9||||||||||||||||||     |   ||||||||||||||||||||||||||||||||||||||
 0123456789012345678901234567890123456789012345678901234567890123
I (339) view: Plot: Length=512, min=-60.000000, max=40.000000
I (339) main: FFT for 1024 complex points take 140472 cycles
I (349) main: End Example.

But I don't know the frequency of the signales. With my application (microphone connected to the ADC) I get:

[2024-07-25_13:32:46:46.542207] W (48745) FFT: Signal audio 64 length
[2024-07-25_13:32:46:46.555478] I (48745) view: Data min[122] = -85.514282, Data max[1] = 11.955421
[2024-07-25_13:32:46:46.556381]  ________________________________________________________________
[2024-07-25_13:32:46:46.557176] 0                                                                |
[2024-07-25_13:32:46:46.570453] 1                                                                |
[2024-07-25_13:32:46:46.571031] 2|                                                               |
[2024-07-25_13:32:46:46.571548] 3                                                                |
[2024-07-25_13:32:46:46.586858] 4                                                                |
[2024-07-25_13:32:46:46.587426] 5 |                                                              |
[2024-07-25_13:32:46:46.603128] 6 |||||                                                          |
[2024-07-25_13:32:46:46.603908] 7      ||||                                                      |
[2024-07-25_13:32:46:46.604542] 8          || |||||||| |  ||                                     |
[2024-07-25_13:32:46:46.618592] 9            |        | ||  ||||||||||||||||||||||||||||||||||||||
[2024-07-25_13:32:46:46.619142]  0123456789012345678901234567890123456789012345678901234567890123
[2024-07-25_13:32:46:46.619679] I (48815) view: Plot: Length=128, min=-60.000000, max=40.000000
[2024-07-25_13:32:46:46.634020] W (48825) FFT: Signal audio 128 length
[2024-07-25_13:32:46:46.634676] I (48825) view: Data min[122] = -85.514282, Data max[1] = 11.955421
[2024-07-25_13:32:46:46.650805]  ________________________________________________________________________________________________________________________________
[2024-07-25_13:32:46:46.667079] 0                                                                                                                                |
[2024-07-25_13:32:46:46.668157] 1                                                                                                                                |
[2024-07-25_13:32:46:46.682961] 2||                                                                                                                              |
[2024-07-25_13:32:46:46.700103] 3                                                                                                                                |
[2024-07-25_13:32:46:46.701265] 4                                                                                                                                |
[2024-07-25_13:32:46:46.714558] 5  |                                                                                                                             |
[2024-07-25_13:32:46:46.730904] 6     |||||||                                                                                                                    |
[2024-07-25_13:32:46:46.746437] 7   ||       || |||||                                                                                                            |
[2024-07-25_13:32:46:46.747402] 8              |     ||||  ||||| || ||| ||   |      ||                                                                           |
[2024-07-25_13:32:46:46.762959] 9                        ||     |  |   |  ||| ||||||  ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
[2024-07-25_13:32:46:46.778493]  01234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567
[2024-07-25_13:32:46:46.779097] I (48975) view: Plot: Length=128, min=-60.000000, max=40.000000
[2024-07-25_13:32:46:46.794639] I (48975) FFT: FFT for 256 complex points take 28091 cycles
[2024-07-25_13:32:46:46.795037] I (48985) FFT: End Example.
[2024-07-25_13:32:49:49.787987] I (51995) FFT: Start Example.

Below is the code I'm using to calculate the FFT of the audio signal:

#include "fft_functions.h"
#include "esp_dsp.h"
#include <math.h>

static const char *FFT_TAG = "FFT";

#define N_SAMPLES 256 //1024
int N = N_SAMPLES;
// Input test array
__attribute__((aligned(16)))
float x1[N_SAMPLES];
__attribute__((aligned(16)))
float x2[N_SAMPLES];
__attribute__((aligned(16)))
float x3[N_SAMPLES];
// Window coefficients
__attribute__((aligned(16)))
float wind[N_SAMPLES];
// working complex array
__attribute__((aligned(16)))
float y_cf[N_SAMPLES * 2];
// Pointers to result arrays
float *y1_cf = &y_cf[0];
float *y2_cf = &y_cf[N_SAMPLES];
float *y3_cf = &y_cf[N_SAMPLES];

// Sum of y1 and y2
__attribute__((aligned(16)))
float sum_y[N_SAMPLES / 2];

void run_fft(float audio_signal[])
{
    esp_err_t fft_ret;
    ESP_LOGI(FFT_TAG, "Start Example.");
    fft_ret = dsps_fft2r_init_fc32(NULL, CONFIG_DSP_MAX_FFT_SIZE);
    if (fft_ret != ESP_OK)
    {
        ESP_LOGE(FFT_TAG, "Not possible to initialize FFT. Error = %i", fft_ret);
        return;
    }
    // Generate hann window
    dsps_wind_hann_f32(wind, N);

    // Convert two input vectors to one complex vector
    for (int i = 0 ; i < N ; i++)
    {
        y_cf[i * 2 + 0] = audio_signal[i] * wind[i];
        y_cf[i * 2 + 1] = 0;
    }
    // FFT
    unsigned int start_b = dsp_get_cpu_cycle_count();
    dsps_fft2r_fc32(y_cf, N);
    unsigned int end_b = dsp_get_cpu_cycle_count();
    // Bit reverse
    dsps_bit_rev_fc32(y_cf, N);
    // Convert one complex vector to two complex vectors
    dsps_cplx2reC_fc32(y_cf, N);

    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);
    }

    ESP_LOGW(FFT_TAG, "Signal audio 64 length");
    dsps_view(y1_cf, N / 2, 64, 10,  -60, 40, '|');
    ESP_LOGW(FFT_TAG, "Signal audio 128 length");
    dsps_view(y1_cf, N / 2, 128/*64*/, 10,  -60, 40, '|');
    ESP_LOGI(FFT_TAG, "FFT for %i complex points take %i cycles", N, end_b - start_b);
    ESP_LOGI(FFT_TAG, "End Example.");
}

In my code run_fft receives as input the array audio_signal. This array contains the ADC values divided by 4095 (12-bit resolution ADC). The values are recovered from the DMA which is set with a sampling frequency of 40 KHz. The DMA is read by an independent task. Below the code I'm using for the ADC:

#include "esp_adc/adc_oneshot.h"
#include "esp_adc/adc_cali.h"
#include "esp_adc/adc_cali_scheme.h"
#include "esp_log.h"
#include "fft_functions.h"
#include "freertos/FreeRTOS.h"
#include "freertos/task.h"
#include <string.h>
//#include <stdint.h>

const static char *TAG = "ADC EXAMPLE";

#include "esp_adc/adc_continuous.h"

#define ADC1_CHAN3          ADC_CHANNEL_3 /* GPIO4 */
#define DEFAULT_ADC_ATTEN       ADC_ATTEN_DB_12 /* ADC_ATTEN_DB_11 */
#define TEST_ADC_CONTINUOUS 1 /* 0 */
#define EXAMPLE_ADC_UNIT                    ADC_UNIT_1
#define _EXAMPLE_ADC_UNIT_STR(unit)         #unit
#define EXAMPLE_ADC_UNIT_STR(unit)          _EXAMPLE_ADC_UNIT_STR(unit)
#define EXAMPLE_ADC_CONV_MODE               ADC_CONV_SINGLE_UNIT_1
#define EXAMPLE_ADC_ATTEN                   ADC_ATTEN_DB_12
#define EXAMPLE_ADC_BIT_WIDTH               SOC_ADC_DIGI_MAX_BITWIDTH

#define EXAMPLE_ADC_OUTPUT_TYPE             ADC_DIGI_OUTPUT_FORMAT_TYPE2
#define EXAMPLE_ADC_GET_CHANNEL(p_data)     ((p_data)->type2.channel)
#define EXAMPLE_ADC_GET_DATA(p_data)        ((p_data)->type2.data)
#define EXAMPLE_READ_LEN                    1024 //256
#define ADC_DMA_ARRAY_SIZE                                  (EXAMPLE_READ_LEN / 4)

float adc_dma_volt_array[ADC_DMA_ARRAY_SIZE] = {0};

static adc_channel_t channel[1] = {ADC_CHANNEL_3};

static TaskHandle_t s_task_handle;
static const char *ADC_TAG = "ADC continuous EXAMPLE";

static bool IRAM_ATTR s_conv_done_cb(adc_continuous_handle_t handle, const adc_continuous_evt_data_t *edata, void *user_data)
{
    BaseType_t mustYield = pdFALSE;

    /* Notify that ADC continuous driver has done enough number of conversions. */
    vTaskNotifyGiveFromISR(s_task_handle, &mustYield);

    return (mustYield == pdTRUE);
}

static void continuous_adc_init(adc_channel_t *channel, uint8_t channel_num, adc_continuous_handle_t *out_handle)
{
    adc_continuous_handle_t handle = NULL;

    adc_continuous_handle_cfg_t adc_config = {
        .max_store_buf_size = 1024,
        .conv_frame_size = EXAMPLE_READ_LEN,
    };
    ESP_ERROR_CHECK(adc_continuous_new_handle(&adc_config, &handle));

    adc_continuous_config_t dig_cfg = {
        .sample_freq_hz = 40000/*20 * 1000*/,
        .conv_mode = EXAMPLE_ADC_CONV_MODE,
        .format = EXAMPLE_ADC_OUTPUT_TYPE,
    };

    adc_digi_pattern_config_t adc_pattern[SOC_ADC_PATT_LEN_MAX] = {0};
    dig_cfg.pattern_num = channel_num;
    for (int i = 0; i < channel_num; i++) {
        adc_pattern[i].atten = EXAMPLE_ADC_ATTEN;
        adc_pattern[i].channel = channel[i] & 0x7;
        adc_pattern[i].unit = EXAMPLE_ADC_UNIT;
        adc_pattern[i].bit_width = EXAMPLE_ADC_BIT_WIDTH;

        ESP_LOGI(ADC_TAG, "adc_pattern[%d].atten is :%"PRIx8, i, adc_pattern[i].atten);
        ESP_LOGI(ADC_TAG, "adc_pattern[%d].channel is :%"PRIx8, i, adc_pattern[i].channel);
        ESP_LOGI(ADC_TAG, "adc_pattern[%d].unit is :%"PRIx8, i, adc_pattern[i].unit);
    }
    dig_cfg.adc_pattern = adc_pattern;
    ESP_ERROR_CHECK(adc_continuous_config(handle, &dig_cfg));

    *out_handle = handle;
}

uint16_t adc1_ch0_value;
double adc_average = 0.0;
adc_oneshot_unit_handle_t adc1_handle;
adc_oneshot_unit_init_cfg_t init_config1 = {
     .unit_id = ADC_UNIT_1,
};
adc_oneshot_chan_cfg_t adc_config = {
     .bitwidth = ADC_BITWIDTH_DEFAULT,
     .atten = DEFAULT_ADC_ATTEN,
};

int adc_raw = 0;
int adc_voltage = 0;

void adc_task(void *pvParameter)
{
    uint32_t ret_num = 0;
    uint8_t result[EXAMPLE_READ_LEN] = {0};
    esp_err_t ret;
    memset(result, 0xcc, EXAMPLE_READ_LEN);

    s_task_handle = xTaskGetCurrentTaskHandle();

    adc_continuous_handle_t handle = NULL;
    continuous_adc_init(channel, sizeof(channel) / sizeof(adc_channel_t), &handle);

    adc_continuous_evt_cbs_t cbs = {
        .on_conv_done = s_conv_done_cb,
    };
    ESP_ERROR_CHECK(adc_continuous_register_event_callbacks(handle, &cbs, NULL));
    ESP_ERROR_CHECK(adc_continuous_start(handle));

    //-------------ADC1 Calibration Init---------------//
    adc_cali_handle_t adc1_cali_handle = NULL;
    bool do_calibration1 = adc_calibration_init(ADC_UNIT_1, DEFAULT_ADC_ATTEN, &adc1_cali_handle);

    while (1)
    {
        /**
            * This is to show you the way to use the ADC continuous mode driver event callback.
            * This `ulTaskNotifyTake` will block when the data processing in the task is fast.
            * However in this example, the data processing (print) is slow, so you barely block here.
            *
            * Without using this event callback (to notify this task), you can still just call
            * `adc_continuous_read()` here in a loop, with/without a certain block timeout.
            */
            ulTaskNotifyTake(pdTRUE, portMAX_DELAY);

            char unit[] = EXAMPLE_ADC_UNIT_STR(EXAMPLE_ADC_UNIT);

            while (1)
            {
                if (current_state != RUNNING) /* For testing this condition can be replaced by 'if (1)' */
                {
                    ret = adc_continuous_read(handle, result, EXAMPLE_READ_LEN, &ret_num, 0);
                    if (ret == ESP_OK)
                    {
                        //ESP_LOGI("TASK", "ret is %x, ret_num is %"PRIu32" bytes", ret, ret_num);
                        //printf("SOC_ADC_DIGI_RESULT_BYTES: %u\n", SOC_ADC_DIGI_RESULT_BYTES);
                        for (int i = 0; i < ret_num; i += SOC_ADC_DIGI_RESULT_BYTES)
                        {
                            adc_digi_output_data_t *p = (adc_digi_output_data_t*)&result[i];
                            uint32_t chan_num = EXAMPLE_ADC_GET_CHANNEL(p);
                            uint32_t data = EXAMPLE_ADC_GET_DATA(p);
                            /* Check the channel number validation, the data is invalid if the channel num exceed the maximum channel */
                            if (chan_num < SOC_ADC_CHANNEL_NUM(EXAMPLE_ADC_UNIT))
                            {
                                //ESP_LOGI(ADC_TAG, "Unit: %s, Channel: %"PRIu32", Value: %"PRIx32, unit, chan_num, data);
                                //printf("%02X\n", data);
                                if (do_calibration1)
                                {
                                    /* Convert the ADC value to voltage and store the value. */
                                    ESP_ERROR_CHECK(adc_cali_raw_to_voltage(adc1_cali_handle, data, &adc_voltage));
                                    adc_dma_volt_array[i / SOC_ADC_DIGI_RESULT_BYTES] = ((float) data) / 4095.0; //(float) adc_voltage;
                                    //ESP_LOGI(ADC_TAG, "ADC%d Channel[%d] Cali Voltage: %d mV", ADC_UNIT_1 + 1, ADC1_CHAN3, adc_voltage);
                                }
                            }
                            else
                            {
                                ESP_LOGW(ADC_TAG, "Invalid data [%s_%"PRIu32"_%"PRIx32"]", unit, chan_num, data);
                            }
                        }

                        //printf("Array ready to run the FFT.\n");
                        run_fft(adc_dma_volt_array);
                        /**
                            * Because printing is slow, so every time you call `ulTaskNotifyTake`, it will immediately return.
                            * To avoid a task watchdog timeout, add a delay here. When you replace the way you process the data,
                            * usually you don't need this delay (as this task will block for a while).
                            */
                        vTaskDelay(3000 /  portTICK_PERIOD_MS/*1*/);
                    }
                    else if (ret == ESP_ERR_TIMEOUT)
                    {
                        //We try to read `EXAMPLE_READ_LEN` until API returns timeout, which means there's no available data
                        break;
                    }
                }
                vTaskDelay(1);
            }
        vTaskDelay(1);
    }

    ESP_ERROR_CHECK(adc_continuous_stop(handle));
    ESP_ERROR_CHECK(adc_continuous_deinit(handle));
}
dmitry1945 commented 1 month ago

Hi @marco711

In dsp_view void dsps_view(const float *data, int32_t len, int width, int height, float min, float max, char view_char); you define amount of points to show on the X axis: width. It means all samples from your input array will be fit to this width points. In the example, it will be N/2, means 128 samples will be fit to the 64 samples.

marco711 commented 1 month ago

Hello @dmitry1945

Thanks for your reply. To be clear, in the following example:

[2024-07-25_13:32:46:46.542207] W (48745) FFT: Signal audio 64 length
[2024-07-25_13:32:46:46.555478] I (48745) view: Data min[122] = -85.514282, Data max[1] = 11.955421
[2024-07-25_13:32:46:46.556381]  ________________________________________________________________
[2024-07-25_13:32:46:46.557176] 0                                                                |
[2024-07-25_13:32:46:46.570453] 1                                                                |
[2024-07-25_13:32:46:46.571031] 2|                                                               |
[2024-07-25_13:32:46:46.571548] 3                                                                |
[2024-07-25_13:32:46:46.586858] 4                                                                |
[2024-07-25_13:32:46:46.587426] 5 |                                                              |
[2024-07-25_13:32:46:46.603128] 6 |||||                                                          |
[2024-07-25_13:32:46:46.603908] 7      ||||                                                      |
[2024-07-25_13:32:46:46.604542] 8          || |||||||| |  ||                                     |
[2024-07-25_13:32:46:46.618592] 9            |        | ||  ||||||||||||||||||||||||||||||||||||||
[2024-07-25_13:32:46:46.619142]  0123456789012345678901234567890123456789012345678901234567890123
[2024-07-25_13:32:46:46.619679] I (48815) view: Plot: Length=128, min=-60.000000, max=40.000000
[2024-07-25_13:32:46:46.634020] W (48825) FFT: Signal audio 128 length
[2024-07-25_13:32:46:46.634676] I (48825) view: Data min[122] = -85.514282, Data max[1] = 11.955421

If the ADC sampling rate is 40 kHz, does this mean that the space between the 64 points in the image is 40 kHz/64 = 0.625 kHz? And if the width is set to 128, will the space between points be 40 kHz/128 = 0.3125 kHz?

dmitry1945 commented 1 month ago

@marco711 If you have sample rate 40k and width = 64, then 40/2/64 = 0.3125 Hz if width is 128, then it's 0.625 Hz

marco711 commented 1 month ago

@dmitry1945 Thank you very much for the explanation.