FFT на STM32F407 Опции

[011119xx]

Перешел на STM32F407, соответственно хочу использовать функции FFT, которые идут с STM32F4xx_DSP_StdPeriph_Lib_V1.0.0. В качестве входных данных использую массив 16-ти битных данных, полученных с АЦП. Судя по всему лучше использовать Real FFT Functions http://tech.munts.com/MCU/Frameworks/AR ... f_f_t.html Но вот каким образом мне не совсем понятно. Смотрел пример, но там вариант для Complex FFT Functions. По подобию примера делаю так:

Код

int16_t analysis_arr[512];      //данные с АЦП
int16_t spec_ana_arr[1024];    

arm_status status;
arm_rfft_instance_q15 S;
arm_cfft_radix4_instance_q15 S_CFFT;
  
status = ARM_MATH_SUCCESS;

/* Initialize the RFFT/RIFFT module */
status = arm_rfft_init_q15(&S, &S_CFFT, 512, 0, 0);

/* Process the data through the RFFT/RIFFT module */
arm_rfft_q15(&S, analysis_arr, spec_ana_arr);

Вопрос к знающим. Правильно ли я делаю?

[scifi]

Видимо да, правильно.
Вот реальный работающий пример:

Код

arm_rfft_instance_f32 rffti;
arm_cfft_radix4_instance_f32 cffti;
arm_rfft_init_f32(&rffti, &cffti, 2048, 0, 1);
arm_rfft_f32(&rffti, f, g);

[011119xx]

Какого размера у вас массивы f и g? Для чего нужно чтобы bitReverseFlag = 1? Можете пояснить почему используется arm_rfft_f32, а не arm_rfft_q31, например?

[scifi]

Какого размера у вас массивы f и g? Для чего нужно чтобы bitReverseFlag = 1? Можете пояснить почему используется arm_rfft_f32, а не arm_rfft_q31, например?

Да забыл:

Код

float32_t f[2048], g[4096];

Про bit reversal немного поискал: по-видимому, на выходе FFT естественным образом получается результат, в котором биты идут в обратном порядке, а с флагом bitReverseFlag = 1 порядок битов снова становится нормальным. Могу и ошибаться, конечно, но работало только так (вывод спектра на график).
Ну а float использовался потому, что так было удобнее при дальнейшей обработке, а времени на вычисления хватало.

[011119xx]

Спасибо.

[011119xx]

Входной массив у вас задан как:

Код

float32_t f[2048]

То есть данные приходящие из АЦП вы предварительно обрабатываете?

[scifi]

То есть данные приходящие из АЦП вы предварительно обрабатываете?

Да.

Код

arm_q15_to_float(adc_data, f, 2048);

[011119xx]

После FFT надо полагать используете обратное преобразование arm_float_to_q15?

[scifi]

После FFT надо полагать используете обратное преобразование arm_float_to_q15?

Нет. Далее так:

Код

arm_cmplx_mag_squared_f32(g, f, 2048);

Это для получения спектра. Потом поиск пиков в спектре:

Код

for (i = 10; i < 500; i++)
{
        if (f[i] > 30.0f)
        {
                ...
        }
}

Ну и так далее.

[011119xx]

Именно так?

Код

arm_cmplx_mag_squared_f32(g, f, 2048);

Судя по примерам там делается так:

Код

arm_cmplx_mag_squared_f32(f, g, 2048);

[scifi]

Уже началось переливание из пустого в порожнее, как мне кажется.
Для справки я использую документацию (в формате HTML), которая шла вместе с CMSIS DSP Library, которая, в свою очередь, установилась при установке ЯРа. Там всё написано.
CMSIS можно скачать отдельно.

[011119xx]

В общем запутался в спектральном анализе. Собственно говоря, делаю так:

Код

float32_t analysis_arr[512];      //данные с АЦП
float32_t spec_ana_arr[1024];
float32_t Output[512];    

arm_status status;
arm_rfft_instance_f32 S;
arm_cfft_radix4_instance_f32 S_CFFT;
  
status = ARM_MATH_SUCCESS;

/* Initialize the RFFT/RIFFT module */
status = arm_rfft_init_f32(&S, &S_CFFT, 512, 0, 1);

/* Process the data through the RFFT/RIFFT module */
arm_rfft_f32(&S, analysis_arr, spec_ana_arr);

arm_cmplx_mag_f32(spec_ana_arr, Output, 512);

Правильно ли я понял, что массив Output содержит амплитуды гармоник? Мне нужно получить амплитуды гармоник с 1-й по 5-ю. Как их найти?

[scifi]

Правильно ли я понял, что массив Output содержит амплитуды гармоник? Мне нужно получить амплитуды гармоник с 1-й по 5-ю. Как их найти?

Похоже, что так.
Ну а гармоники находятся на известных частотах? Если так, то пересчитываем частоты в индексы массива Output[] и смотрим соответствующие элементы.
Для отладки необходимо выводить на график весь спектр, чтобы убедиться, что он выглядит как ожидается, и отловить все глупые ошибки.

[011119xx]

Ну а гармоники находятся на известных частотах? Если так, то пересчитываем частоты в индексы массива Output[] и смотрим соответствующие элементы.

Думаю, что гармоники на известных частотах, но как пересчитать частоты в индексы массива Output[]? И что делать если окажется, что гармоники на неизвестных частотах?

[scifi]

Думаю, что гармоники на известных частотах, но как пересчитать частоты в индексы массива Output[]?

Почитайте про дискретное преобразование Фурье. Там шаг по частоте - это 1/T или что-то такое, где T - общее время оцифровки.

И что делать если окажется, что гармоники на неизвестных частотах?

Искать нужную частоту. Например, путём поиска в спектре компонент, превышающих заданный порог.

[011119xx]

Порассуждаем. T - общее время оцифровки. В моем случае 512 отсчетов, между отсчетами 5 мкс, итого Т = (512 - 1) * 5 мкс = 2,555 мс. f = 1/T = 391 Гц. То есть это частота первой гармоники f1. Частота второй гармоники f2 = 2* f1 = 782 Гц. И так далее. Теперь смотрим в массив Output[], который представляет собой спектр амплитуд. Он состоит из 512 элементов. Чтобы вычислить амплитуду первой гармоники Am1 нужно выполнить следующее:

Код

am1 = 0;
for(j = 0; j < 512; j++)                                    
    am1 = am1 + Output[j] * sin(3.14 * j) / (512 - 1));

Am1 = (2 * am1 / (512 - 1));

Для второй гармоники:

Код

am2 = 0;
for(j = 0; j < 512; j++)                                    
    am2 = am2 + Output[j] * sin(2 * 3.14 * j) / (512 - 1));

Am2 = (2 * am2 / (512 - 1));

Обобщенно для первых пяти гармоник:

Код

float32_t Am[5];

for(i = 0; i < 5; i++)
{
    am = 0;
    for(j = 0; j < 512; j++)                                    
        am = am + Output[j] * sin(i * 3.14 * j) / (512 - 1));

    Am[i] = (2 * am / (512 - 1));
}

Это верно?

[011119xx]

Понадобилось провести линейную интерполяцию над данными, полученными из АЦП. В библиотеке есть пример. Пытаюсь понять как же его использовать. В начале идет инициализация структуры:

Код

arm_linear_interp_instance_f32 S = {188495, -3.141592653589793238, 0.00005f, &arm_linear_interep_table[0]};

Структура определена как:

Код

  typedef struct
  {
    uint32_t nValues;
    float32_t x1;
    float32_t xSpacing;
    float32_t *pYData;          /**< pointer to the table of Y values */
  } arm_linear_interp_instance_f32;

nValues - размер массива данных, *pYData - адрес массива данных, xSpacing - это видимо шаг по оси Х, в случае массива - шаг индекса массива, x1 - вот это не понятно что такое.
Вопрос в том, как правильно инициализировать эту структуру, иначе дальнейшие действия не имеют смысла?

[011119xx]

Вновь возвращаюсь с вопросом. Необходимо провести FFT над данными из АЦП. Часть кода такова:

Код

int16_t analysis_arr[512];      //данные с АЦП
int16_t spec_ana_arr[1024];
int16_t amplitude_arr[512];    

arm_status status;
arm_rfft_instance_q15 S;
arm_cfft_radix4_instance_q15 S_CFFT;
  
status = ARM_MATH_SUCCESS;

status = arm_rfft_init_q15(&S, &S_CFFT, 512, 0, 1);

arm_rfft_q15(&S, analysis_arr, spec_ana_arr);

arm_cmplx_mag_q15(spec_ana_arr, amplitude_arr, 512);

В функции arm_cmplx_mag_q15 выполняется:

Код

    /* C[0] = sqrt(A[0] * A[0] + A[1] * A[1]) */
    real = *pSrc++;
    imag = *pSrc++;
    acc0 = __SMUAD(real, real);
    acc1 = __SMUAD(imag, imag);
    /* store the result in 2.14 format in the destination buffer. */
    arm_sqrt_q15((q15_t) (((q63_t) acc0 + acc1) >> 17), pDst++);

Вопрос в том, для чего нужно сдвигать сумму квадратов аж на 17 разрядов?

Еще вопросы. После выполнения arm_rfft_q15 исходный массив analysis_arr меняется. Почему? Выходной формат функции arm_rfft_q15 9.7, что это означает?

[nowhow]

Вновь возвращаюсь с вопросом. Необходимо провести FFT над данными из АЦП. Часть кода такова:

real = *pSrc++;
imag = *pSrc++;
acc0 = __SMUAD(real, real);
acc1 = __SMUAD(imag, imag);
/* store the result in 2.14 format in the destination buffer. */
arm_sqrt_q15((q15_t) (((q63_t) acc0 + acc1) >> 17), pDst++);[/code]
Вопрос в том, для чего нужно сдвигать сумму квадратов аж на 17 разрядов?

Еще вопросы. После выполнения arm_rfft_q15 исходный массив analysis_arr меняется. Почему? Выходной формат функции arm_rfft_q15 9.7, что это означает?

сдвиг нужен для приведения форматов чисел.
9.7 означает, что 9 бит отводится под целую часть числа, а 7 бит под дробную.
посмотрите описание функций, там должны быть входные и выходные форматы данных.

[Sekat]

Видимо да, правильно.
Вот реальный работающий пример:

Код

arm_rfft_instance_f32 rffti;
arm_cfft_radix4_instance_f32 cffti;
arm_rfft_init_f32(&rffti, &cffti, 2048, 0, 1);
arm_rfft_f32(&rffti, f, g);

Запускаю так:

Код

arm_q15_to_float(testInput_q15_512, f, 512);
arm_rfft_instance_f32 rffti;
arm_cfft_radix4_instance_f32 cffti;
arm_rfft_init_f32(&rffti, &cffti, 2048, 0, 1);
arm_rfft_f32(&rffti, f, g);
arm_cmplx_mag_squared_f32(g, f, 512);

То есть делаю FFT для real f32, как привел пример scifi - все великолепно. Вижу хорошую палку на 10 индексе. Тестовые отсчеты сгенерированы в MatLab в формате q15. 512 - отсчетов, 10 периодов синуса, целочисленно укладывающихся в заданное окно. Размах от 0xC000 до 0x3FFF

Тоже самое пытаюсь сделать для real q15

Код

arm_rfft_instance_q15 rffti;
arm_cfft_radix4_instance_q15 cffti;
arm_rfft_init_q15( &rffti, &cffti, 512, 0, 1);
arm_rfft_q15( &MyS, testInput_q15_512, testOutput);
arm_cmplx_mag_q15(testOutput, amplitude_arr, 512);

- полный бред на выходе, больше похожий на белый шум.
Если задать третью строчку в коде так: arm_rfft_init_q15( &rffti, &cffti, 512, 1, 1); - то есть обратное преобразование
вместо прямого, то картина становится лучше. Появляется не очень внятный и размазанный пик на нужном месте http://electronix.ru/forum/style_emoticons...t/smile3046.gif .
Тестовые отсчеты сгенерированы в MatLab в формате q15. 512 - отсчетов, 10 периодов синуса, целочисленно укладывающихся в заданное окно.

Голова сломана, выходные проходят бездарно! дошел до перекомпиляции последней версии CMSIS DSP. Результат тот же!
Направьте на истинный путь !

[Allregia]

Чтобы не плодить тем, спрошу тут.
Я только начал разбираться с DSP функциями в F407, поэтому не ругайте больно за детские вопросы.
1) У меня входные данные - 24 бита с АЦП. Это какому формату в библиотеке соответствует?
2) 31 тар FIR - как прикинуть макс. время выполнения в F407 на штатной тактовой?
3) безотносительно FIR - как на этом проце максимально быстро преобразовать переменную типа char задом наперед ? (т.е. реверс битов сделать).
В других процах я байт разбивал на 2 тетрады, и каждую тетраду по таблице (чуть медленнее чем с одной таблицей, но зато таблица не 256 а всего 16 байт), у F407 памяти много, можно и не разбивать.

[KnightIgor]

3) безотносительно FIR - как на этом проце максимально быстро преобразовать переменную типа char задом наперед ? (т.е. реверс битов сделать).

Можно изобразить на ассемблере, используя инструкцию RBIT. Она меняет местами биты в 32-битном регистре, но очевидно, что после операции самый старший байт регистра будет содержать реверсированный младший. Объявление функции может выглядеть как:

uint8_t reverse_8bits(uint8_t value);

а код в ассемблерном файле (с учетом соглашения на передачу параметров, когда R0 - есть первый параметр):

Код

reverse_8bits  PROC
              EXPORT reverse_8bits  [WEAK]

              RBIT  R0, R0
              REV   R0, R0
              BX    LR

              ENDP

P.S. Код на работоспособность сам еще не проверял.
P.P.S. Проверил под симулятором: работает.

Сообщение отредактировал KnightIgor - Mar 2 2013, 13:07

[HHIMERA]

реверс битов сделать

Код

b = __REV(__RBIT(a));

[Allregia]

Код

b = __REV(__RBIT(a));

Мерси!!! Работает изумительно.

P.S. мне это надо было для управления дисплеем по SPI через USART (USART в отличие от SPI умеет только один порядок битов. а дисплею надо наоборот).

[hd44780]

Ламерский вопрос - какую либу подключать -
arm_cortexM4lf_math.lib (Little endian and Floating Point Unit on Cortex-M4)
или
arm_cortexM4bf_math.lib (Big endian and Floating Point Unit on Cortex-M4) ?

Никак не врублюсь, какой индеец в этом STM32F407

Или есть у кого нормальный пример, как из выборки АЦП спектр получить?

Спасибо.

[_Артём_]

Никак не врублюсь, какой индеец в этом STM32F407

Little...по умолчанию...

[hd44780]

Little...по умолчанию...

Спасибо.
Подключил пока arm_cortexM4lf_math.lib, посмотрю, что будет ...

По мотивам прочитанного здесь насочинял такое (но ещё не закочил):

CODE

// выходной спектр FFT
uint8_t SpectrumOut [N_WAVE_HALF/2];

// N_WAVE_HALF = 2048
static q15_t t [ N_WAVE_HALF ];
static float32_t f [ N_WAVE_HALF ], g [ N_WAVE_HALF * 2 ];

// Выполнение FFT
// pInData - входные данные (выборка FFT) - N_WAVE_HALF
void DoFFT ( uint8_t *pInData )
{
arm_rfft_instance_f32 rffti;
arm_cfft_radix4_instance_f32 cffti;
int i;

// 1. Подготовка входных данных
for ( i = 0; i < N_WAVE_HALF; i ++ )
t [ i ] = pInData [i];

// 1-й арг - int32_t
arm_q15_to_float ( t, f, N_WAVE_HALF );

// 2. Инициализация FFT
arm_rfft_init_f32 ( &rffti, &cffti, N_WAVE_HALF, 0, 1 );

// 3. FFT
arm_rfft_f32 ( &rffti, f, g );

// 4. Спектр
arm_cmplx_mag_squared_f32(g, f, 2048);
} // DoFFT

Будет ли работать - не знаю ....

[hd44780]

Получилось вот такое:

CODE

bool DoFFT ( uint8_t *pInData )
{
arm_rfft_instance_f32 rffti;
arm_cfft_radix4_instance_f32 cffti;
arm_status stat;
int i;

// 1. Подготовка входных данных
for ( i = 0; i < N_WAVE; i ++ )
t [ i ] = pInData [i];

// 1-й арг - int32_t t -> f
arm_q15_to_float ( t, f, N_WAVE );

// 2. Инициализация FFT
stat = arm_rfft_init_f32 ( &rffti, &cffti, N_WAVE, 0, 1 );

if ( stat != ARM_MATH_SUCCESS )
return false;

// 3. FFT f -> g
arm_rfft_f32 ( &rffti, f, g );

// 4. Спектр g -> f
arm_cmplx_mag_squared_f32 ( g, f, N_WAVE );

for ( i= 0; i < N_WAVE/2; i ++)
{
SpectrumOut[i] = (uint8_t) f[i];
} // fof

return true;
} // DoFFT

В результате в массиве f почти все числа меньше 1, соответственно в выходном целочисленном SpectrumOut одни нули
Помогите, кто может.
Входной массив, поступающий в функцию нормальный. По нему осциллограмму рисую, огна нормальная.

Спасибо.

Сообщение отредактировал IgorKossak - Mar 6 2013, 07:27

Причина редактирования: [codebox] для длинного кода, [code] - для короткого!!!