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[])

[prig]

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

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

Вероятно, Вы пыталисьсь обобщить частный случай реализации целочисленного БПФ и сделали ошибку.

...
Тут теперь уже Ваши объяснения очень тяжело понять, так что даже дискутировать не могу. Единственно результат сдвигать надо вправо в последнем предложении, а округлить надо до сдвига.
К стати: вопрос на засыпку для разминки мозгов: почему при умножении двух 16-разрядных знаковых чисел в допкоде результат получается 32-разрядным, когда по казалось бы очевидной логике должен быть 31 разряд? )))
...

- Речь была о реализации функции умножения 1.15*1.15 -> 1.15 с помощью целочисленного умножения со знаком (ADI ещё используют описание типа 16.0*16.0 -> 32.0):

"Результат операции целочисленного умножения со знаком аргументов в формате 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 и потом округлить."

А вот как говорят в ADI о реализации этой функции для одной из инструкций Блэкфина:
Signed fraction. Multiply 1.15 * 1.15 to produce 1.31 results after left-shift correction.
Round 1.31 format value at bit 16. (RND_MOD bit in the ASTAT register controls the rounding.)
Saturate the result to 1.15 precision in destination register half.
Result is between minimum -1 and maximum 1-2-15 (or, expressed in hex, between minimum 0x8000 and maximum 0x7FFF).

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

Возможно, Вы пытались поведать о чём-то другом, подразумевая ваше понимание какого-то частного случая.
Аналогично, заданный Вами вопрос о разрядности тоже оказался не слишком корректен.
Для такой детерминированной функции, как целочисленное умножение со знаком (16.0*16.0 -> 32.0), можно говорить о диапазоне значений результата или его порядке, а также о формате результата и диапазоне значений этого формата. И в этой части всё очевидно.

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

- Вероятно, Вы моделировали с табличным синусом или чем-то простеньким. Попробуйте на реальном сигнале типа OFDM. Фиксированная постадийная нормализация (может осуществляется в процессе вычисления бабочки) приведёт к потере гейна и, соответственно, к увеличению шума самого преобразования. И в отдельных случаях потери будут чувствительные (например, для 16-бит 8К БПФ на сигнале DVB-T это уже хорошо заметно, проверялось на реальном же железе, а о 16-ти битах для 32К DVB-T2 даже говорить нечего).

Так что, ваше обобщение касается только частных случаев. В общем случае оно неверно.



П.С. Судя по всему, Вы всё время пытаетесь интерпретировать единственный частный случай реализации целочисленного БПФ.
Это не есть правильно. Таки, лучше подучить матчасть для общего случая.

[Krys]

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

- Подчёркнутое утверждение неверно. Вращающиеся вектора - неправильное объяснение.

И что же неправильно во вращающихся векторах?

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

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

Вероятно, Вы пыталисьсь обобщить частный случай реализации целочисленного БПФ и сделали ошибку.

Я привёл условия типа "если каждую квадратуру выхода БПФ считать пронормированным к 1". Я согласен, что это для частного случая. Но для него я не сделал ошибку.

Т.е., два старших бита результата умножения 1.15*1.15 в формате 2.30 будут всегда знаковыми. Т.о., для получения операции 1.15*1.15 -> 1.15 можно арифметически сдвинуть влево результат 2.30 и потом округлить."
А вот как говорят в ADI о реализации этой функции для одной из инструкций Блэкфина:
Signed fraction. Multiply 1.15 * 1.15 to produce 1.31 results after left-shift correction.
Round 1.31 format value at bit 16. (RND_MOD bit in the ASTAT register controls the rounding.)
Saturate the result to 1.15 precision in destination register half.
Result is between minimum -1 and maximum 1-2-15 (or, expressed in hex, between minimum 0x8000 and maximum 0x7FFF).

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

Согласен, если рассматривать иностранную цитату непогрешимой, то моё утверждение неверно. Но если считать, что цитата не является неоспоримой аксиомой, то ошибся не я, а цитата и Вы вместе с ней. Жирным с подчёркиванием я выделил места под вопросом.

"Всегда" - ошибка. Именно про это был мой "вопрос на засыпку", и засыпка произошла.

"Saturate the result to 1.15" - в рамках предыдущего "всегда" сатурация не требуется, там просто нечего сатурировать.

"after left-shift correction" - мне непонятно, почему влево. Вот пример:
Есть переменная 8 битов с форматом 4.4, её нужно вместить с потерей дробной части в 4-битную переменную 4.0. Что мы делаем:

Код

исходная переменная, физическая разрядность (не так, которую мы
подразумеваем форматом 4.4):
                            7 6 5 4 3 2 1 0

Надо вместить в переменную:         3 2 1 0

Если просто скопировать биты 3..0 в биты 3..0 - сохранится младшая часть

Тогда исходную переменную надо сдвинуть вправо:

                            0 0 0 0 7 6 5 4
                            
А затем просто скопировать биты 3..0 в биты 3..0:
                            
Исходная переменная:        0 0 0 0 7 6 5 4
Копируем:                           | | | |
Конечная переменная:                7 6 5 4

Сдвиг-то вправо.

Аналогично, заданный Вами вопрос о разрядности тоже оказался не слишком корректен.

Ещё как корректен, просто это нужно прощупать "на своей шкуре".

Для такой детерминированной функции, как целочисленное умножение со знаком (16.0*16.0 -> 32.0), можно говорить о диапазоне значений результата или его порядке, а также о формате результата и диапазоне значений этого формата. И в этой части всё очевидно.

Очевидно, но первое, что очевидно, оказывается неверным. Вопрос на засыпку был "почему при умножении двух 16-разрядных знаковых чисел в допкоде результат получается 32-разрядным, когда по казалось бы очевидной логике должен быть 31 разряд"? Очевидно это 31 разряд, а не 32. Потому что при умножении полезные данные сосредоточены в младших 15 битах операндов, следовательно 15 * 15 получаем 30 и 1 бит знака. 31 разряд в результате. Откуда 32?

- Вероятно, Вы моделировали с табличным синусом или чем-то простеньким.

Согласен, тут я сделал свой вывод только для своего сигнала. Он был реальный, но с сильным преобладанием несущей, на фоне которой нужно было обнаружить отражения с допплером.

П.С. Судя по всему, Вы всё время пытаетесь интерпретировать единственный частный случай реализации целочисленного БПФ.
Это не есть правильно. Таки, лучше подучить матчасть для общего случая.

А Вы уверены, что смогли охватить этот самый общий случай? Он бесконечен. Все случаи учесть невозможно. Так-то хорошо все утверждения парировать, что это не годится, в общем случае работать не обязано. Универсальный козырь.