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]

Мы говорим про функцию, которая рассчитывает амплитуду частотных составляющих по их мнимой и действительной частям?


Ну я Вам про форматы пытался выше объяснить. Если что - спрашивайте конкретные вопросы, поясню более детально.
В рамках конкретной функции такие форматы можно объяснить: если каждую квадратуру выхода БПФ считать пронормированным к 1, то амплитуда в корень из 2 будет больше 1, поэтому при сохранении входного формата возникнет переполнение, следовательно положение точки нужно сдвинуть.
Но на самом деле этого не требуется конкретно для БПФ, т.к. заранее известно, что комплексные числа - это вращающиеся вектора, которые не могут выходить за единичную окружность, тогда и модуль не может превышать 1.


А надо-то вправо. БПФ усилил на гейн. Надо вернуть обратно.


А частота дискретизации тогда какая?


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

Да, arm_cmplx_mag_q15() рассчитывает амплитуду по мнимой и комплексной части.

По поводу форматов конкретный вопрос собственно такой: У меня есть знаковое целочисленное число - отсчеты моего АЦП. Нужно ли его преобразовать, чтобы получить формат 1.15 ибо именно в таком формате мне нужно подавать данные на вход FFT? И вот на выходе arm_cmplx_mag_q15() формат 2.14, как мне его преобразовать, чтобы получить значение амплитуды в отсчетах АЦП?

По поводу гейна и сдвига. В том-то и дело, если сдвигать ВПРАВО, то в выходном массиве после вызова arm_cmplx_mag_q15() одни нули. Если не сдвигать, то картинка характеристики "обрезается", нет "нижних" маленьких "тычков" (гармоники с малой амплитудой).

Частота дискретизации 1024 Гц. Повторил эксперимент для частоты сгенерированного математически сигнала 32 Гц (т.е. 1024 кратно 32), амплитуда 600. Картинка исходного сигнала вот:



На выходе вот такая картинка:

[Krys]

По поводу форматов конкретный вопрос собственно такой: У меня есть знаковое целочисленное число - отсчеты моего АЦП. Нужно ли его преобразовать, чтобы получить формат 1.15 ибо именно в таком формате мне нужно подавать данные на вход FFT? И вот на выходе arm_cmplx_mag_q15() формат 2.14, как мне его преобразовать, чтобы получить значение амплитуды в отсчетах АЦП?

Положение точки в числах с фиксированной точкой - понятие чисто условное, которое существует только в головах разработчиков, чтобы понимать, какое реальное число представлено в коде из скажем 8 битов.
Например, дана переменная 8 битов со значением битов 48. Сказано, что эта переменная имеет формат s3.5. Какое же реальное число закодировано в этих 8 битах? Закодировано 1,5. А если это s4.4 - закодировано 3.
Внутри вычислителя положение точки никак не фигурирует, для него это всегда полностью целые числа, без дробной части. Но разработчик, зная, что входной 8-битный знаковый поток он представил в уме, скажем, как нормализованный к 1 (т.е. s1.7 - это реальные числа сдвинутые влево на 7 битов), понимает, что результат ему надо задвинуть обратно на те же 7 битов, чтобы получить правильный результат в его реальных числах, которые он представлял в уме.
Т.е. у нас есть инструмент, который оперирует только целыми числами. И у нас есть реальные числа, которые никак не укладываются в диапазон этого инструмента. Нам надо сначала в уме растянуть наши числа на полную шкалу инструмента, а затем результат (он опять в целых числах) обратно сжать на столько же битов. Для 8-битных чисел, например, инструмент понимает числа до 127, а реальные числа у нас не больше 1. Поэтому все наши реальные числа надо умножить на 127, чтобы влезть в полную шкалу инструмента.
Возвращаясь к ответу на исходный вопрос:
На входе предполагается (в уме) формат 1.15, а на выходе 2.14. Как это понять?
Отвлечённый пример: пусть на входе формат 1.15 и на выходе формат 1.15. Это означает, что положение точки в одном и том же месте. Если на входе мы свои числа растянули на полную шкалу, сдвинув влево на 15 разрядов, то на выходе нужно сдвинуть вправо. А если мы на входе ничего не сдвигали (это всё равно мысленная операция в конце концов), то и на выходе ничего сдвигать не надо.
Теперь вернёмся к конкретно нашим цифрам: на входе 1.15, на выходе 2.14. Даже если мы вход мысленно не сдвигали, то выход будет всё равно сдвинут функцией на 1 бит, при том вправо. Чтобы получить число в том же масштабе, что и входное, потребуется сдвинут результат на 1 бит влево.
Вот Вам и ответ. Т.е. вход сдвигать не надо вообще, а выход - на 1 бит влево.

По поводу гейна и сдвига. В том-то и дело, если сдвигать ВПРАВО, то в выходном массиве после вызова arm_cmplx_mag_q15() одни нули. Если не сдвигать, то картинка характеристики "обрезается", нет "нижних" маленьких "тычков" (гармоники с малой амплитудой).

Понятно, Вы методом научного тыка решили, что сдвигать всё же надо )))
Я придерживаюсь мнения, что сдвигать вообще не надо никуда, функция БПФ и так всё сделала за Вас, оно и логично по-идее с точки зрения разработки юзер-френдли интерфейса, чтобы пользователь не парил голову, сколько стадий бабочки прошло, ему нужен только готовый результат.
Однако вот тут мне непонятно:

Если не сдвигать, то картинка характеристики "обрезается", нет "нижних" маленьких "тычков" (гармоники с малой амплитудой).

Откуда у Вас после сдвига влево спектр приобретает "нижние маленькие тычки", если без сдвига их нет? До сдвига эта информация должна быть уже заложена в несдвинутый сигнал. А если она заложена, то как она обрезается отсутствием каких-либо действий (случай без сдвига)? Не должно такого быть. Подумайте, может, где-то ошибка.

Частота дискретизации 1024 Гц. Повторил эксперимент для частоты сгенерированного математически сигнала 32 Гц (т.е. 1024 кратно 32), амплитуда 600. Картинка исходного сигнала вот:
На выходе вот такая картинка:

Другое дело. Вот с этого сразу и надо было начинать.
К стати, по картинке спектра: левая и правая палки имеют разный уровень. Этого быть не должно. Уж не знаю, откуда взялось. Может, функция кривая, может, мнимая часть где-то пролезла (хотя почти не реально). Может, другие участники форума подскажут, откуда берётся этот эффект. Если не найдётся причина - предлагаю входные числа затолкать в матлаб один к одному, и над ними проделать fft(). Там должно быть одинаково.

[prig]

...
Теперь вернёмся к конкретно нашим цифрам: на входе 1.15, на выходе 2.14. Даже если мы вход мысленно не сдвигали, то выход будет всё равно сдвинут функцией на 1 бит, при том вправо. Чтобы получить число в том же масштабе, что и входное, потребуется сдвинут результат на 1 бит влево.
Вот Вам и ответ. Т.е. вход сдвигать не надо вообще, а выход - на 1 бит влево.
...

Звиздец какое объяснение.
Я так понимаю, Вы пытаетесь донести вариант реализации умножения 1.15*1.15 -> 1.15 с использованием целочисленного умножителя.

Однако, к вопросу "И вот на выходе arm_cmplx_mag_q15() формат 2.14, как мне его преобразовать, чтобы получить значение амплитуды в отсчетах АЦП?" это имеет только косвенное отношение. И вот почему:

Диапазон значений корня из суммы квадратов вещественной и мнимой части в формате 1.15 (диапазон каждой что-то около -0.99997...0.99997) будет выходить за границы диапазона -1...1. Т.е., в общем случае, для представления модуля комплексного числа в вышеуказанном формате 16-ти битным числом потребуется фомат 2.14 (примерно -1.99994...1.99994). С округлением, ессно. Как именно получается этот результат, вопрос десятый.
С учётом того, что ТС собирается использовать нормировочный коэффициент (см. вопрос), использование формата 2.14 как 1.15 это всего лишь вопрос значения этого нормировочного коэффициента. Результат этой функции никуда сдвигать не надо.

[Krys]

Звиздец какое объяснение.

Осторожнее с эмоциями, здесь им не место )

Я так понимаю, Вы пытаетесь донести вариант реализации умножения 1.15*1.15 -> 1.15 с использованием целочисленного умножителя.

Не обязательно умножения. К стати, с умножением всё сложнее. Я скорее пытался объяснить, что значит такой формат в виде двух чисел с точкой, и почему он всего лишь мысленный.

Диапазон значений корня из суммы квадратов вещественной и мнимой части в формате 1.15 (диапазон каждой что-то около -0.99997...0.99997) будет выходить за границы диапазона -1...1. Т.е., в общем случае, для представления модуля комплексного числа в вышеуказанном формате 16-ти битным числом потребуется фомат 2.14 (примерно -1.99994...1.99994). С округлением, ессно. Как именно получается этот результат, вопрос десятый.

Я там выше подобный механизм расписывал для бабочки:

В рамках конкретной функции такие форматы можно объяснить: если каждую квадратуру выхода БПФ считать пронормированным к 1, то амплитуда в корень из 2 будет больше 1, поэтому при сохранении входного формата возникнет переполнение, следовательно положение точки нужно сдвинуть.
Но на самом деле этого не требуется конкретно для БПФ, т.к. заранее известно, что комплексные числа - это вращающиеся вектора, которые не могут выходить за единичную окружность, тогда и модуль не может превышать 1.

С учётом того, что ТС собирается использовать нормировочный коэффициент (см. вопрос), использование формата 2.14 как 1.15 это всего лишь вопрос значения этого нормировочного коэффициента. Результат этой функции никуда сдвигать не надо.

А такое объяснение не понял теперь я )))

[prig]

Осторожнее с эмоциями, здесь им не место )



Не обязательно умножения. К стати, с умножением всё сложнее. Я скорее пытался объяснить, что значит такой формат в виде двух чисел с точкой, и почему он всего лишь мысленный.


Я там выше подобный механизм расписывал для бабочки:



А такое объяснение не понял теперь я )))

Осторожнее с изобретением объяснений типа "на пальцах". Можно попасть впросак.
Арифметика с фикс. точкой не так банальна, как это может показаться, это да.
Но все эти вопросы давно и очень хорошо изложены в соответствующей литературе.
Соответственно, доморощенные объяснения могут как прийти в противоречие с устоявшейся терминологией и практикой, так и привести к ошибкам.

С умножением всё достаточно просто.
Достаточно осознать, что из себя представляют диапазоны соответствующих форматов с фикс.точкой и диапазоны результатов операций.
Вот несколько наиболее показательных примеров:
Порядок диапазона для беззнаковых аргументов 16-бит и результата целочисленного умножения без знака (32 бита): 2^16 и 2^32. Как бы очивидно.
То же самое, но для аргументов со знаком: 2^16 и 2^31.
Теперь для 1.15.
Результат операции целочисленного умножения со знаком аргументов в формате 1.15 имеет тот же порядок диапазона - 2^31. Диапазон значений результата укладывается в -1...1.
Формат результата - 2.30. Порядок диапазона формата 2.30 - 2^32. Диапазон значений формата 2.30 укладывается в -2...2.
Т.е., два старших бита результата умножения 1.15*1.15 в формате 2.30 будут всегда знаковыми. Т.о., для получения операции 1.15*1.15 -> 1.15 можно арифметически сдвинуть влево результат 2.30 и потом округлить.

С БПФ у Вас как раз незадача и вышла.
Если использовать арифметику 1.15*1.15 -> 1.15 и 1.15+-1.15 -> 1.15 , возможность переполнения при вычислении БПФ гарантирована.
Именно из-за этого используется 1 или 2 "защитных" бита для радикса-2 или радикса-4 соответстенно.
Т.е. весь входной массив нормируется до необходимых значений. Нормировочный коэффициент запоминается.
Кроме того, на таком формате обычно используется групповая плавающая точка. Т.е., постадийная нормировка по-необходимости. Нормировочный к-т обновляется, ессно.
С промежуточной арифметикой 2.30+-2.30 -> 2.30 всё немного по-другому, но постадийная нормировка нормировка всё равно обычно используется.

1.15 и 2.14 - это фактически одно и то же. Различать их позволяет как раз нормировочный коэффициент.