Форум разработчиков электроники ELECTRONIX.ru

Помощь - Поиск - Пользователи - Календарь

Полная версия этой страницы: FFT на STM32F407

Форум разработчиков электроники ELECTRONIX.ru > Микроконтроллеры (MCs) > ARM

011119xx

Apr 17 2012, 05:15

Перешел на 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

Apr 17 2012, 06:00

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

Код

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

Apr 17 2012, 06:23

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

scifi

Apr 17 2012, 06:44

Цитата(011119xx @ Apr 17 2012, 10:23)

Да забыл:

Код

float32_t f[2048], g[4096];

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

011119xx

Apr 17 2012, 06:51

Спасибо.

011119xx

Apr 17 2012, 08:37

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

Код

float32_t f[2048]

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

scifi

Apr 17 2012, 08:50

Цитата(011119xx @ Apr 17 2012, 12:37)

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

Да.

Код

arm_q15_to_float(adc_data, f, 2048);

011119xx

Apr 17 2012, 09:31

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

scifi

Apr 17 2012, 10:33

Цитата(011119xx @ Apr 17 2012, 13:31)

После 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

Apr 17 2012, 11:11

Именно так?

Код

arm_cmplx_mag_squared_f32(g, f, 2048);

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

Код

arm_cmplx_mag_squared_f32(f, g, 2048);

scifi

Apr 17 2012, 11:23

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

011119xx

Apr 24 2012, 08:11

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

Код

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

Apr 24 2012, 09:12

Цитата(011119xx @ Apr 24 2012, 12:11)

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

011119xx

Apr 24 2012, 09:43

Цитата(scifi @ Apr 24 2012, 14:12)

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

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

scifi

Apr 24 2012, 09:59

Цитата(011119xx @ Apr 24 2012, 13:43)

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

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

Цитата(011119xx @ Apr 24 2012, 13:43)

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

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

011119xx

Apr 25 2012, 08:34

Порассуждаем. 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

Apr 26 2012, 04:42

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

Код

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

Jun 14 2012, 04:51

Вновь возвращаюсь с вопросом. Необходимо провести 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

Jul 4 2012, 06:56

Цитата(011119xx @ Jun 14 2012, 08:51)

Вновь возвращаюсь с вопросом. Необходимо провести 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

Mar 2 2013, 09:46

Цитата(scifi @ Apr 17 2012, 10:00)

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

Код

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

Mar 2 2013, 10:42

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

KnightIgor

Mar 2 2013, 12:54

Цитата(Allregia @ Mar 2 2013, 11:42)

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. Проверил под симулятором: работает.

HHIMERA

Mar 2 2013, 14:18

Цитата(Allregia @ Mar 2 2013, 13:42)

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

Код

b = __REV(__RBIT(a));

Allregia

Mar 3 2013, 12:09

Цитата(HHIMERA @ Mar 2 2013, 16:18)

Код

b = __REV(__RBIT(a));

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

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

hd44780

Mar 3 2013, 18:36

Ламерский вопрос - какую либу подключать -
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

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

Спасибо.

_Артём_

Mar 3 2013, 20:52

Цитата(hd44780 @ Mar 3 2013, 20:36)

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

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

hd44780

Mar 4 2013, 07:31

Цитата(_Артём_ @ Mar 3 2013, 22:52)

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

Mar 4 2013, 19:24

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

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 одни нули

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

Спасибо.

Для просмотра полной версии этой страницы, пожалуйста, пройдите по ссылке.