FFT stm32f407 лишние гармоники в выходном массиве Опции

[yanvasiij]

Доброго времени суток!

Процессор stm32f407, библиотека DSP от Stm32CubeMx. Делаю выборку из 1024 отсчетов встроенного АЦП (12 бит) с частотой 1024 Гц. Когда набирается 1024 отсчета отправляю данные в функции быстрого преобразования Фурье из либы DSP.

Код в более приятном виде можно посмотреть ЗДЕСЬ.

CODE

arm_cfft_radix2_instance_q15 fftInstance;
static u32 curInputCursor=0;
static u32 dataIsReady = 0; /**< Флаг для начала преобразования фурье */
static q15_t input [2048]; /**< входной массив данных */
static q15_t output [1024]; /**< выходной массив данных */

/** @brief прерывание по таймеру с частотой 1024 Гц */
extern "C" void timerInterrupt (void)
{
input[curInputCursor++]=uhADCxConvertedValue; /**< Набираю 1024 замера */
if (curInputCursor>=1024)
{
curInputCursor = 0;
dataIsReady = 1; /**< Флаг, сигнализирующий, что выборка готова */
HAL_TIM_Base_DeInit(&htim1); /**< Как только набралось отключаю таймер, чтобы остановить забор замеров */
}
invertLed(); /**< Инвертирую ногу, чтобы на осциллографе проверить действительно ли частота дискретизации 1024 Гц (да, именно так) */
}

/** @brief задача, которая занимается преобразованием фурье */
static void fftTask (void *p)
{
q15_t maxValue; /* Max FFT value is stored here */
uint32_t maxIndex; /* Index in Output array where max value is */

while(1)
{
if (dataIsReady)
{
dataIsReady = 0;

arm_cfft_radix2_init_q15 (&fftInstance, 1024, 0, 1);
arm_cfft_radix2_q15(&fftInstance, input);
arm_cmplx_mag_q15 (input, output, 1024);
u32 index;
saveFft();
timInit(); /**< Снова запускаю таймер для подготовки следующих 1024 замеров */
}
vTaskDelay(5000);
}
}

Если я все правильно понимаю то в массиве output[] должна получиться частотная характеристика с шагом 1 Гц. При этом output[0] - это постоянная составляющая, все остальное гармоники.
Подаю синусоидальный сигнал 50 Гц с амплитудой 1 вольт, смещенный на 1 вольт. Ожидаю, что в output[0] будет значение постоянной составляющей, а в output[50] амплитуда 50-ой гармоники.
Однако в результате output[0] всегда равен нулю, в output[1] некоторое значение очень похожее на постоянную составлющую, в output[3] тоже непонятно какое значение (181), а в output[50] как и ожидалось значение гармоники (оно смещается если изменить частоту синусоидального сигнала). Вдобавок нет нет да и появятся непонятные значения в разных гармониках (то output[200], то output[511] и т.п. совершенно случайным образом). Добавил RC-цепочку на вход - не помогло.

Буду очень признателен, если кто растолкует, что я делаю неправильно.

Вот здесь картинка с характеристикой (строю на графике массив output[])

[yanvasiij]

Одной из причин может быть то, что подаваемая частота 50Гц не строго синхронна частоте выборок АЦП, поэтому возникает размытие спектра. Попробуйте тест вашего БПФ не на реальных данных с АЦП, а на заранее подготовленном в Matlab буфере данных, где частота сигнала строго синхронна получится. Можно даже в симуляторе прогнать наверное.

Вообщем начал проверять, как Вы и сказали (сгенерировал массив с синусоидой и беру данные из него). Получилось тоже самое. Но обнаружил ошибку: в массиве input[] данные нужно размещать вот так: input[n] = действительная часть, input[n+1]=комплексная часть (n=0,2,4,8,10...1022). В моем случае, если я правильно понимаю, комплексная часть входного массива должна быть равна нулю. Вдобавок, как выяснилось из документации на либу, данные во входном массиве к функции arm_cmplx_mag_q15() нужно располагать в каком-то особом формате (вот тут описание функции, в официальной доке). Вот, что значит "The function implements 1.15 by 1.15 multiplications and finally output is converted into 2.14 format", убей не пойму? Полазил по форумам, нашел вот тут, что надо предварительно смещать на 10 бит влево. Сместил, совершенно не понимая зачем. В итоге код принял вот такой вид:

Вот тут с подсветкой

CODE

#define FFT_LEN 1024

u32 uhADCxConvertedValue;
static u32 curInputCursor=0;
static u32 dataIsReady = 0;

arm_cfft_radix2_instance_q15 fftInstance;
static q15_t input [FFT_LEN*2];
static q15_t output [FFT_LEN];

/**
* @brief Обработчик прерывания от таймера, тактирующего измерения (частота 1024 Гц)
*
* @return
*/
extern "C" void timerInterrupt (void)
{
input[curInputCursor++]=uhADCxConvertedValue; /**< действительная часть */
input[curInputCursor++] = 0; /**< Комплексная часть */

/** Если набралось FFT_LEN измерений, то останавливаю таймер */
if (curInputCursor>=FFT_LEN*2)
{
curInputCursor = 0;
dataIsReady = 1;
HAL_TIM_Base_DeInit(&htim1);
}
invertLed(); /**< инвертирую ногу, чтобы проверить частоту дискретизации */
}

/**
* @brief Запустить FFT c сохранением данных на SD-карте
*
* @param frequency
* @param fileName
*/
static void runFft (TimerFrequencyType frequency, const char *fileName)
{
timInit(frequency); /**< Инициалзирую таймер с частой 1024 Гц */
while (dataIsReady == 0) taskYIELD(); /**< Жду готовности данных */
dataIsReady = 0;

/** Собственно CFFT */
arm_cfft_radix2_init_q15 (&fftInstance, FFT_LEN, 0, 1);
arm_cfft_radix2_q15(&fftInstance, input);
/** Смещению на 10 бит */
for (u32 i=0; i<2048; i++) input[i] = input[i]<<10;
/** Вычисляю амплитуды гармоник */
arm_cmplx_mag_q15 (input, output, FFT_LEN);
/** Вывожу в файл и в веб */
saveFft(frequency, fileName);
}

Получилось вот такая картинка (на реальном сигнале с АЦП):



Мне так и не понятно правильно ли я понял, что нужно смещать данные во входном массиве на 10 бит влево прежде чем отправлять в arm_cmplx_mag_q15()? Если да, то почему??? И во вторых, в каком формате данные в выходном массиве output[]? Это амплитуды гармоник в отсчетах АЦП или в дБ? Разъясните, пожалуйста.

- Для "фирменных" библиотек опыт показывает, что если правильные данные данные в правильном формате правильно размещены во входном массиве, результат соответствует теории с точностью до эффектов насыщения и т.д. Намёк ясен?
- О спектре оконной функции у Вас есть представление? Если нет, то Вы явно не с того начали.

- А я и не выставлял претензий к "фирменным" библиотекам и в их работоспособности не сомневаюсь. Вопрос вообщем в том и состоит, правильно ли я размещаю данные и в том ли формате? Потому и привел код и попытался подробно объяснить, что я делаю, может я где ошибся или чего неправильно понял. Вопрос именно в этом и был.
- Если речь идет спектре сигнала ограниченного во времени, и вообще о том как должен выглядеть спектр на выходе с ДПФ, то да у меня есть об этом представление.