Библиотеки DSP - RFFT - странные результаты

Недавно я пытался выполнить вычисления БПФ на моей оценочной плате STM32F4-Discovery, а затем отправить их на ПК. Я изучил свою проблему - мне кажется, что я делаю что-то не так с функциями БПФ, предоставленными производителем.

Я использую библиотеки CMSIS-DSP. На данный момент я генерировал образцы с кодом (если это работает правильно, я сделаю выборку с помощью микрофона).

Я использую arm_rfft_fast_f32, поскольку в будущем мои данные будут плавающими, но результаты, которые я получаю в моем выходном массиве, безумны (я думаю) - я получаю частоты ниже 0.

number_of_samples = 512; (l_probek in code)
dt = 1/freq/number_of_samples

Вот мой код

float32_t buffer_input[l_probek];
uint16_t i;
uint8_t mode;
float32_t dt;
float32_t freq;
bool DoFlag = false;
bool UBFlag = false;
uint32_t rozmiar = 4*l_probek;

union
{
    float32_t f[l_probek];
    uint8_t b[4*l_probek];
}data_out;


union
{
    float32_t f[l_probek];
    uint8_t b[4*l_probek];
}data_mag;

union
{
    float32_t f;
    uint8_t b[4];
}czest_rozdz;


/* Pointers ------------------------------------------------------------------*/
arm_rfft_fast_instance_f32 S;
arm_cfft_radix4_instance_f32 S_CFFT;
uint16_t output;
/* ---------------------------------------------------------------------------*/
int main(void)
{
    freq = 5000;
    dt = 0.000000390625;


    _GPIO();
    _LED();
    _NVIC();    
    _EXTI(0);

    arm_rfft_fast_init_f32(&S, l_probek);
    GPIO_SetBits(GPIOD, LED_Green);

    mode = 2;


    //----------------- Infinite loop
  while (1)
    {
        if(true)//(UBFlag == true)

                    for(i=0; i<l_probek; ++i)
                    {
                        buffer_input[i] = (float32_t) 15*sin(2*PI*freq*i*dt);
                    }

            //Obliczanie FFT
            arm_rfft_fast_f32(&S, buffer_input, data_out.f, 0);
            //Obliczanie modulow
            arm_cmplx_mag_f32(data_out.f, data_mag.f, l_probek);

            USART_putdata(USART1, data_out.b, data_mag.b, rozmiar);
            //USART_putdata(USART1, czest_rozdz.b, data_mag.b, rozmiar);
            GPIO_ToggleBits(GPIOD, LED_Orange);
            //mode++;
            //UBFlag = false;

        }

    }
}

c arm cortex-m fft
person Jejh    schedule 17.02.2017    source источник
comment
Вы проверили, что ваши входные образцы подходят для вашего теста? Кроме того, каково значение l_probek? Это 512?   -  person Dave S    schedule 17.02.2017
comment
@DaveS Похоже, он рассчитывает тестовую синусоиду в buffer_input.   -  person tofro    schedule 17.02.2017
comment
По поводу этой линии - как вы определили амплитуду 15? buffer_input [i] = (float32_t) 15 * sin (2 * PI freq i * dt);   -  person Dave S    schedule 17.02.2017

Ответы (2)


arrow_upward
7
arrow_downward

Я использую arm_rfft_fast_f32, поскольку в будущем мои данные будут плавающими, но результаты, которые я получаю в моем выходном массиве, безумны (я думаю) - я получаю частоты ниже 0.

Функция arm_rfft_fast_f32 не возвращает частоты, но довольно сложные коэффициенты, вычисленные с использованием быстрого преобразования Фурье (БПФ). Таким образом, вполне разумно, чтобы эти коэффициенты были отрицательными. В частности, ожидаемые коэффициенты для входного однопериодного sin тестового сигнала с амплитудой 15 будут следующими:

0.0,     0.0; // special case packing real-valued X[0] and X[N/2]
0.0, -3840.0; // X[1]
0.0,     0.0; // X[2]
0.0,     0.0; // X[3]
...
0.0,     0.0; // X[255]

Обратите внимание, что, как указано в документации, первая два выхода соответствуют чисто реальным коэффициентам X[0] и X[N/2] (вы должны быть особенно осторожны с этим особым случаем при последующем вызове arm_cmplx_mag_f32; см. последний пункт ниже).

Частота каждого из этих частотных компонентов задается k*fs/N, где N - это количество выборок (в вашем случае l_probek), а fs = 1/dt - частота выборки (в вашем случае freq*l_probek):

X[0] -> 0*freq*l_probek/l_probek =              0
X[1] -> 1*freq*l_probek/l_probek =   freq =  5000
X[2] -> 2*freq*l_probek/l_probek = 2*freq = 10000
X[3] -> 3*freq*l_probek/l_probek = 2*freq = 15000
...

Наконец, из-за особой упаковки первых двух значений вам нужно быть осторожным при вычислении N/2+1 звездных величин:

// General case for the magnitudes
arm_cmplx_mag_f32(data_out.f+2, data_mag.f+1, l_probek/2 - 1);
// Handle special cases
data_mag.f[0]          = data_out.f[0];
data_mag.f[l_probek/2] = data_out.f[1];
person SleuthEye    schedule 17.02.2017
comment
Большое спасибо SleuthEye. Блин, я начинаю думать, что я безмозглый;) Еще я напортачил с генерацией семплов для моего БПФ. Неправильно принял dt на выбранную мной частоту. Также сегодня утром я проверил, что было бы разумнее использовать функцию sinf для генерации значений. - person Jejh; 18.02.2017
comment
У меня к вам еще один вопрос - нормально ли, что после вызова функции arm_rfft_fast_f3 мой входной массив изменяется? - person Jejh; 20.02.2017
comment
Да, вкратце, как работает эта функция. Более длинный ответ заключается в том, что arm_rfft_fast_f32 использует a>, который выполняет вычисления на месте с использованием входного буфера (таким образом изменяя входные данные). - person SleuthEye; 20.02.2017
comment
Еще одна вещь - CORTEX-M4 (например, stm32F4) uC с прямым порядком байтов или прямым порядком байтов по умолчанию? - person Jejh; 20.02.2017
comment
Я считаю, что руководство пользователя упоминает прямой порядок байтов доступа к памяти. См. Также этот ответ. - person SleuthEye; 20.02.2017

arrow_upward
1
arrow_downward

В продолжение вышеупомянутого ответа, что потрясающе, некоторые дополнительные пояснения, на которые мне потребовалось время, чтобы понять.

  • Элементы разрешения по частоте центрируются на целевой частоте, поэтому, например, в приведенном выше примере X [0] представляет от -2500 Гц до 2500 Гц с центром в нуле, X [1] составляет от 2500 Гц до 7500 Гц с центром в 5000 Гц и т. Д.

  • Обычно интерполяция частот внутри интервала осуществляется путем просмотра энергии соседних интервалов (см. https://dspguru.com/dsp/howtos/how-to-interpolate-fft-peak/), если вы сделаете это, вам нужно будет убедиться, что ваш массив величин достаточно велик для бины + Найквиста и что бин над Найквистом равен 0, но обратите внимание, что многие методы интерполяции требуют комплексных значений (например, Куинн, Якобсон), поэтому убедитесь, что вы интерполировали, прежде чем находить величины.

  • Приведенный выше код особого случая работает, потому что нет комплексного компонента значений DC и Найквиста, и, следовательно, величина - это просто действительная часть.

  • Однако в приведенном выше коде есть ошибка - хотя мнимые части компонентов DC и Найквиста всегда равны нулю, реальная часть все еще может быть отрицательной, поэтому вам нужно взять абсолютное значение, чтобы получить величину:

// Handle special cases
data_mag.f[0]          = fabs(data_out.f[0]);
data_mag.f[l_probek/2] = fabs(data_out.f[1]);
person Andy Piper    schedule 09.08.2019