Помощь - Поиск - Пользователи - Календарь
Полная версия этой страницы: FFT stm32f407 лишние гармоники в выходном массиве
Форум разработчиков электроники ELECTRONIX.ru > Цифровая обработка сигналов - ЦОС (DSP) > Алгоритмы ЦОС (DSP)
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[])

Нажмите для просмотра прикрепленного файла
Krys
Одной из причин может быть то, что подаваемая частота 50Гц не строго синхронна частоте выборок АЦП, поэтому возникает размытие спектра. Попробуйте тест вашего БПФ не на реальных данных с АЦП, а на заранее подготовленном в Matlab буфере данных, где частота сигнала строго синхронна получится. Можно даже в симуляторе прогнать наверное.
prig
Цитата(yanvasiij @ Aug 19 2015, 13:48) *
...
Однако в результате output[0] всегда равен нулю, в output[1] некоторое значение ...
...
Вдобавок нет нет да и появятся непонятные значения в разных гармониках (то output[200], то output[511] и т.п. совершенно случайным образом).
...

- Для "фирменных" библиотек опыт показывает, что если правильные данные данные в правильном формате правильно размещены во входном массиве, результат соответствует теории с точностью до эффектов насыщения и т.д. Намёк ясен?
- О спектре оконной функции у Вас есть представление? Если нет, то Вы явно не с того начали.
Krys
Не умничайте )) Напишите, что конкретно
Alex11
Конкретно - если не наложить окно на измеренный сигнал перед Фурье, то можно узнать о себе много нового и интересного в зависимости от погоды в Лондоне и др. внешних факторов.
Lmx2315
Цитата(Alex11 @ Aug 21 2015, 14:01) *
Конкретно - если не наложить окно на измеренный сигнал перед Фурье, то можно узнать о себе много нового и интересного в зависимости от погоды в Лондоне и др. внешних факторов.

..по умолчанию окно всегда наложено - прямоугольное, что как бы ничего не объясняет.
з.ы.
я бы сделал 1024 отчёта табличного синуса и посмотрел как он выглядит, только синус и всё, потом он и постоянка (искусственная const).
Потом два синуса, потом синус на рабочей частоте.
Krys
Ну я что и предлагал - тесты на заранее подготовленных данных, а не на "живых" ))
yanvasiij
Цитата(Krys @ Aug 21 2015, 11:23) *
Одной из причин может быть то, что подаваемая частота 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[]? Это амплитуды гармоник в отсчетах АЦП или в дБ? Разъясните, пожалуйста.

Цитата(prig @ Aug 21 2015, 12:26) *
- Для "фирменных" библиотек опыт показывает, что если правильные данные данные в правильном формате правильно размещены во входном массиве, результат соответствует теории с точностью до эффектов насыщения и т.д. Намёк ясен?
- О спектре оконной функции у Вас есть представление? Если нет, то Вы явно не с того начали.


- А я и не выставлял претензий к "фирменным" библиотекам и в их работоспособности не сомневаюсь. Вопрос вообщем в том и состоит, правильно ли я размещаю данные и в том ли формате? Потому и привел код и попытался подробно объяснить, что я делаю, может я где ошибся или чего неправильно понял. Вопрос именно в этом и был.
- Если речь идет спектре сигнала ограниченного во времени, и вообще о том как должен выглядеть спектр на выходе с ДПФ, то да у меня есть об этом представление.
Alex11
Может быть мне следовало изъясняться тщательнЕе. Естественно, прямоугольное окно есть всегда, только оно не обеспечивает отсутствие разрывов функции и ее производной на краях окна. В результате, при периодичности сигнала не кратной размеру окна, получаем кучу артефактов. На живом сигнале это бывает всегда в той или иной степени.
yanvasiij
Цитата(Alex11 @ Aug 21 2015, 16:44) *
Может быть мне следовало изъясняться тщательнЕе. Естественно, прямоугольное окно есть всегда, только оно не обеспечивает отсутствие разрывов функции и ее производной на краях окна. В результате, при периодичности сигнала не кратной размеру окна, получаем кучу артефактов. На живом сигнале это бывает всегда в той или иной степени.


Я, если честно не нашел ни одного упоминая о том, использует ли DSP stm32f4xx библиотека оконное сглаживание. Выходит, что если нет, то подобные артефакты могут появится? Если да, то как это можно проверить? Сгенерировать сигнал не кратный размеру окна?
Lmx2315
Цитата(yanvasiij @ Aug 21 2015, 15:43) *
Мне так и не понятно правильно ли я понял, что нужно смещать данные во входном массиве на 10 бит влево прежде чем отправлять в arm_cmplx_mag_q15()? Если да, то почему??? И во вторых, в каком формате данные в выходном массиве output[]? Это амплитуды гармоник в отсчетах АЦП или в дБ? Разъясните, пожалуйста.

..10 влево - это 1024 раза, что подозрительно коррелирует с размером вашей БПФ.
думаю в выходном массиве точно не Дб, а линейные величины.
yanvasiij
Цитата(Lmx2315 @ Aug 21 2015, 16:57) *
..10 влево - это 1024 раза, что подозрительно коррелирует с размером вашей БПФ.
думаю в выходном массиве точно не Дб, а линейные величины.


Да, вот только значение в output[50] (тот пик, что на картинке) равно 10543, а у меня АЦП 12 бит, т.е. максимум 4095. Как так?
prig
Цитата(Lmx2315 @ Aug 21 2015, 13:54) *
..по умолчанию окно всегда наложено - прямоугольное, что как бы ничего не объясняет.
з.ы.
я бы сделал 1024 отчёта табличного синуса и посмотрел как он выглядит, только синус и всё, потом он и постоянка (искусственная const).
Потом два синуса, потом синус на рабочей частоте.


- Если накладыватель не в курсе, считайте, что не наложено (последствия описаны Alex11).
- Если у ТС "выбивает" постоянку, что ему дадут табличные синусы?
Скорее всего, ТС что-то не то или не туда грузит на входе или не там смотрит по выходу (о чём я как бы совсем уж прозрачно "намекал").
М.б. битреверс там отдельной функцией, и он не приложился им в нужном месте.

Крче, ТС следовало бы сперва доки по библиотеке как следует вкурить. Ну и окна ...

Цитата(yanvasiij @ Aug 21 2015, 14:52) *
Я, если честно не нашел ни одного упоминая о том, использует ли DSP stm32f4xx библиотека оконное сглаживание. Выходит, что если нет, то подобные артефакты могут появится? Если да, то как это можно проверить? Сгенерировать сигнал не кратный размеру окна?

И не найдёте. Оконная функция - это очень простая операция (или отсутствие оной в случае прямоугольного окна).
Выбор оконной функции далеко не элементарен и полностью лежит на совести разработчика.
Lmx2315
Цитата(yanvasiij @ Aug 21 2015, 16:04) *
Да, вот только значение в output[50] (тот пик, что на картинке) равно 10543, а у меня АЦП 12 бит, т.е. максимум 4095. Как так?

..ну и что? интерес практический или академический?
ваши 4095 были во временной области, а вы уже смотрите частотную , где совсем другие размерности.
Если вам надо сравнивать сигналы между собой - отнормируйте спектр по известному сигналу.
Подайте виртуальный синус размером с дин. диапазон вашего АЦП и посмотрите сколько попугаев он будет занимать в частотной области.
Кстати, тут уже будут сильно сказываться применяемые вами окна.
Alex11
Кроме всего изложенного выше, в спектре для определения амплитуды исходного сигнала нельзя использовать одиночный отсчет. Нужно брать квадратный корень из суммы квадратов отсчетов, попадающих в полосу сигнала. Ширина этой полосы определяется оконной функцией и требуемой точностью.
Krys
Цитата(Lmx2315 @ Aug 21 2015, 19:49) *
Если вам надо сравнивать сигналы между собой - отнормируйте спектр по известному сигналу.
Подайте виртуальный синус размером с дин. диапазон вашего АЦП и посмотрите сколько попугаев он будет занимать в частотной области.

А математически вычислить никак? )) По крайней мере не для живого сигнала и для кратных частот. Я думаю, с этого надо начать.
А для живого сигнала, конечно, уже, как написано тут:
Цитата(Alex11 @ Aug 21 2015, 23:03) *
Кроме всего изложенного выше, в спектре для определения амплитуды исходного сигнала нельзя использовать одиночный отсчет. Нужно брать квадратный корень из суммы квадратов отсчетов, попадающих в полосу сигнала. Ширина этой полосы определяется оконной функцией и требуемой точностью.
yanvasiij
Цитата(Lmx2315 @ Aug 21 2015, 17:49) *
..ну и что? интерес практический или академический?
ваши 4095 были во временной области, а вы уже смотрите частотную , где совсем другие размерности.
Если вам надо сравнивать сигналы между собой - отнормируйте спектр по известному сигналу.
Подайте виртуальный синус размером с дин. диапазон вашего АЦП и посмотрите сколько попугаев он будет занимать в частотной области.
Кстати, тут уже будут сильно сказываться применяемые вами окна.



Цитата(Alex11 @ Aug 21 2015, 21:03) *
Кроме всего изложенного выше, в спектре для определения амплитуды исходного сигнала нельзя использовать одиночный отсчет. Нужно брать квадратный корень из суммы квадратов отсчетов, попадающих в полосу сигнала. Ширина этой полосы определяется оконной функцией и требуемой точностью.



Цитата(Krys @ Aug 24 2015, 07:15) *
А математически вычислить никак? )) По крайней мере не для живого сигнала и для кратных частот. Я думаю, с этого надо начать.
А для живого сигнала, конечно, уже, как написано тут:


Правильно ли я понимаю что речь идет просто о вычислении амплитуды комплексного числа?
Нажмите для просмотра прикрепленного файла
Но ведь именно это и делает arm_cmplx_mag_q15(). После ДПФ input[] уже содержит в себе действительную и мнимую составляющие, из них можно вычислить и амплитуду и фазу. Я вычислил амплитуды и построил по ним график, но результат мне не понятен. Или я что-то путаю?
Krys
Цитата(yanvasiij @ Aug 21 2015, 18:43) *
Вообщем начал проверять, как Вы и сказали (сгенерировал массив с синусоидой и беру данные из него).
А какую частоту синусоиды взяли? Можете файл прикрепить с отсчётами? Фактически интересует, сколько сэмплов период повторения. Но файл лучше )



Цитата(yanvasiij @ Aug 21 2015, 18:43) *
Но обнаружил ошибку: в массиве input[] данные нужно размещать вот так: input[n] = действительная часть, input[n+1]=комплексная часть (n=0,2,4,8,10...1022). В моем случае, если я правильно понимаю, комплексная часть входного массива должна быть равна нулю.
А в чём Ваша ошибка была? Вы вообще неправильно данные на вход подали? Да, для вещественного синуса комплексная часть равна нулю.



Цитата(yanvasiij @ Aug 21 2015, 18:43) *
Вот, что значит "The function implements 1.15 by 1.15 multiplications and finally output is converted into 2.14 format", убей не пойму?
Почитайте на вики про бабочку в операциях БПФ, тогда все вопросы исчезнут. В 2х словах: формат 1.15 означает 16-битное представление, 1 бит под целую часть (он же знак, поэтому под целую часть нет ничего, и число должно быть меньше 1), а 15 битов - под дробную часть. Грубо говоря, в таком формате все числа пронормированы к 1.
Бабочка выполняет умножение входного числа (отсчёта синуса) на поворачивающий множитель. Поворачивающий множитель по определению не превышает 1, поэтому его формат тоже выбран 1.15 (хотя, это чуть-чуть некорректно, и в своей реализации я выбрал честный формат). После умножения бабочка делает сложение. Отсюда разрядность должна увеличиться на 1 бит. Разрядная сетка у нас фиксирована 16 битами, поэтому приходится отбросить младший бит дробной части, чтобы не словить переполнение. Отсюда получается результирующий формат на выходе бабочки 2.14. Бабочка выполняется над массивом данных 1024 сэмпла 10 раз, поэтому в результате точка съезжает аж на 10 разрядов вправо, про это говорят "гейн БПФ".


Цитата(yanvasiij @ Aug 21 2015, 18:43) *
Вот, что значит "The function implements 1.15 by 1.15 multiplications and finally output is converted into 2.14 format", убей не пойму? Полазил по форумам, нашел вот тут, что надо предварительно смещать на 10 бит влево. Сместил, совершенно не понимая зачем.

Чтобы получить правильную амплитуду, нужно результат поделить на вышеупомянутый гейн. Можно и не делить физически, а "держать в уме" наличие усиленного результата. Так что делить или не делить - от этого качественно картинка меняться не должна, кривая спектра всё равно будет выглядеть одинаково внешне, и палки частот всё равно должны быть в нужных местах. Вы сначала добейтесь правильного положения палок, а затем уже думайте о правильном измерении абсолютного уровня амплитуды палки.
yanvasiij
Цитата(Krys @ Aug 24 2015, 09:00) *
А какую частоту синусоиды взяли? Можете файл прикрепить с отсчётами? Фактически интересует, сколько сэмплов период повторения. Но файл лучше )


А в чём Ваша ошибка была? Вы вообще неправильно данные на вход подали? Да, для вещественного синуса комплексная часть равна нулю.



Почитайте на вики про бабочку в операциях БПФ, тогда все вопросы исчезнут. В 2х словах: формат 1.15 означает 16-битное представление, 1 бит под целую часть (он же знак, поэтому под целую часть нет ничего, и число должно быть меньше 1), а 15 битов - под дробную часть. Грубо говоря, в таком формате все числа пронормированы к 1.
Бабочка выполняет умножение входного числа (отсчёта синуса) на поворачивающий множитель. Поворачивающий множитель по определению не превышает 1, поэтому его формат тоже выбран 1.15 (хотя, это чуть-чуть некорректно, и в своей реализации я выбрал честный формат). После умножения бабочка делает сложение. Отсюда разрядность должна увеличиться на 1 бит. Разрядная сетка у нас фиксирована 16 битами, поэтому приходится отбросить младший бит дробной части, чтобы не словить переполнение. Отсюда получается результирующий формат на выходе бабочки 2.14. Бабочка выполняется над массивом данных 1024 сэмпла 10 раз, поэтому в результате точка съезжает аж на 10 разрядов вправо, про это говорят "гейн БПФ".



Чтобы получить правильную амплитуду, нужно результат поделить на вышеупомянутый гейн. Можно и не делить физически, а "держать в уме" наличие усиленного результата. Так что делить или не делить - от этого качественно картинка меняться не должна, кривая спектра всё равно будет выглядеть одинаково внешне, и палки частот всё равно должны быть в нужных местах. Вы сначала добейтесь правильного положения палок, а затем уже думайте о правильном измерении абсолютного уровня амплитуды палки.



Чуть позже смогу выложить файл с синусоидой. Ошибка была в том, что я ложил данные во входной массив "подряд", а нужно было писать ноль в ячейки с комплексными значениями.
Я добился того, что "палки" находятся "на своих местах". Про бабочку почитаю подробно, пока только и знаю что весь принцип БПФ основан на бабочках (они могут быть по основанию 2,4,8, теоретически возможно и более, но практически слабо реализуемо). А что касается гейна БПФ - ведь я и сместил на 10 разрядов (или если хотите разделил на 2^10), но и после этого на выходе значения даже не близкие к реальным значениям амплитуд гармоник. Выход тоже пронормирован к 1 верно? Значит, чтобы получить амплитуды реальные амплитуды гармоник, то мне нужно сделать так: output[i]/2^15 * 2^12 (мой АЦП 12 бит)?
Krys
Цитата(yanvasiij @ Aug 24 2015, 10:39) *
Правильно ли я понимаю что речь идет просто о вычислении амплитуды комплексного числа?
Я уже сам запутался, кому и что надо )) Вы напишите ещё раз, что хотите рассчитать и что не срастается?


Цитата(yanvasiij @ Aug 24 2015, 10:39) *
Но ведь именно это и делает arm_cmplx_mag_q15(). После ДПФ input[] уже содержит в себе действительную и мнимую составляющие, из них можно вычислить и амплитуду и фазу. Я вычислил амплитуды и построил по ним график, но результат мне не понятен. Или я что-то путаю?

А вообще, Вам бы книжку почитать по БПФ. Мне нравится пример на стр. 151 в книжке Лайонс_-_ЦОС_2е_издание_русский_перевод_2006 (можно скачать в местных закромах). Там хотя бы понятно, как образуются амплитуды и отсчёты, откуда берутся такие числа. Хорошо всё разжёвано. У Вас должно получиться также, как в примере. Да и вообще в этой книжке хорошо даётся материал касаемо данного вопроса.
yanvasiij
Цитата(Krys @ Aug 24 2015, 11:35) *
Я уже сам запутался, кому и что надо )) Вы напишите ещё раз, что хотите рассчитать и что не срастается?



А вообще, Вам бы книжку почитать по БПФ. Мне нравится пример на стр. 151 в книжке Лайонс_-_ЦОС_2е_издание_русский_перевод_2006 (можно скачать в местных закромах). Там хотя бы понятно, как образуются амплитуды и отсчёты, откуда берутся такие числа. Хорошо всё разжёвано. У Вас должно получиться также, как в примере. Да и вообще в этой книжке хорошо даётся материал касаемо данного вопроса.



- Что хочу получить: амплитуды гармоник в отсчетах АЦП.

- И да, пожалуй, пойду перечитаю про БПФ. Спасибо за источник!
Krys
Цитата(yanvasiij @ Aug 24 2015, 11:56) *
А что касается гейна БПФ - ведь я и сместил на 10 разрядов (или если хотите разделил на 2^10)
Я не программист, поэтому в Ваши исходники копаться не полезу, спрашиваю так: вход с АЦП у Вас как переменная описан с плавающей точкой или с фиксированной? Каким конкретно типом, какой разрядности? И что значит "сместил на 10 разрядов", какой операцией? Сместил вход или выход БПФ?

Цитата(yanvasiij @ Aug 24 2015, 11:56) *
но и после этого на выходе значения даже не близкие к реальным значениям амплитуд гармоник.
Покажите, пожалуйста, полный получившийся спектр. Вы приводили выше картинку, но она похоже неполная, там нет отрицательной частоты. И вот Вы пишете, что "значения даже не близкие". Напишите прямо: подал столько-то, поэтому по расчётам (как рассчитывал?) ожидаю получить столько-то, но по рисунку спектра видно столько-то.

Цитата(yanvasiij @ Aug 24 2015, 11:56) *
Выход тоже пронормирован к 1 верно? Значит, чтобы получить амплитуды реальные амплитуды гармоник, то мне нужно сделать так: output[i]/2^15 * 2^12 (мой АЦП 12 бит)?
А Вы на вход подавали просто 12-битные отсчёты АЦП? Тогда если результат поделить на 1024, то должна получиться половинная амплитуда входной синусоиды в попугаях входного 12-битного сигнала.
yanvasiij
Цитата(Krys @ Aug 24 2015, 12:41) *
Я не программист, поэтому в Ваши исходники копаться не полезу, спрашиваю так: вход с АЦП у Вас как переменная описан с плавающей точкой или с фиксированной? Каким конкретно типом, какой разрядности? И что значит "сместил на 10 разрядов", какой операцией? Сместил вход или выход БПФ?

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

А Вы на вход подавали просто 12-битные отсчёты АЦП? Тогда если результат поделить на 1024, то должна получиться половинная амплитуда входной синусоиды в попугаях входного 12-битного сигнала.


Вход АЦП - целочисленная перменная (12 бит), которая хранится в uint16. Эти значения накапливаются в массиве input[] (он имеет тип q15_t, который по сути просто переопределение int16) без каких бы то нибыл изменений (просто взял АЦП и переложил в этот массив) в таком порядке: input[n] - действительная часть, input[n+1] - комплексная (n=0,2,4,6,8...) - в таком виде их нужно там располагать, если я правильно понял документацию на DSP библиотеку. Вот в действительную часть я положил значения АЦП. Этот массив отправляется в функцию arm_cfft_radix2_q15(), которая вообщем то и делает комплексное быстрое преобразование Фурье по основанию 2. Результат в комплексном виде кладет во все тот же массив input[]. Далее, чтобы вычислить амплитуду комплексного числа, в библиотеке DSP есть функция arm_cmplx_mag_q15(), которая в выходном массиве для каждого комплексного числа (в мое случае пары членов input[n] и input[n+1]) находит значение амплитуды (корень квадратный из суммы действительного и мнимого компонента). Однако, чтобы эту функция правильно посчитала, нужно, как я понял, во входном массиве input[] каждый член сдвинуть влево на 10 бит: for (u32 i=0; i<2048; i++) input[i] = input[i]<<10 и только потом передавать функции arm_cmplx_mag_q15() в качестве параметра. В результате я рассчитывал получить массив амплитуд гармоник, составляющий мой сигнал. Вообщем-то массив и получился, с "палкой" в том месте, где я и рассчитывал получить его (output[0] и output[50] - ведь я подал смещенную синусоиду с частотой 50 Гц). Если Вы спрашивали значения этого выходного массива, то вот они в десятичной системе счисления. Подал синусоиду 500 милливольт, при разрядности АЦП 12 бит я рассчитывал увидеть число 610 в десятичной системе (или 262h хексе), а получил 17438 в десятичной (или 441E в хексе). В полном виде на картинке он выглядит вот так:

Нажмите для просмотра прикрепленного файла
Krys
Спасибо, теперь всё гораздо понятнее.
А входной сигнал выложите, пожалуйста.

И Вы пишете, что подали 500мВ, а какое напряжение полной шкалы? АЦП двухполярный? Т.е. 1 бит под знак, а значения в доп.коде на выходе? И вообще, лучше давайте продолжать работать не с живыми данными АЦП. Сначала надо отладить всё на идеальной модели входного сигнала, где всё точно.

С какой разрядностью на входе и выходе работает функция arm_cmplx_mag_q15()?

И ещё Вы пишете:
Цитата
Однако, чтобы эту функция правильно посчитала, нужно, как я понял, во входном массиве input[] каждый член сдвинуть влево на 10 бит
Если не будет выхода за предела разрядной сетки, сдвиг можно сделать хоть до, хоть после функции. Но Вы ещё не назвали разрядность, поэтому как будет известно - будем думать.

И сдвиг вправо, а не влево. Это ошибка в Ваших расчётах или просто опечатка? Если ошибка, то прошу все вышеприведённые цифры повторить заново после исправления.
yanvasiij
Цитата(Krys @ Aug 24 2015, 14:34) *
Спасибо, теперь всё гораздо понятнее.
А входной сигнал выложите, пожалуйста.

И Вы пишете, что подали 500мВ, а какое напряжение полной шкалы? АЦП двухполярный? Т.е. 1 бит под знак, а значения в доп.коде на выходе? И вообще, лучше давайте продолжать работать не с живыми данными АЦП. Сначала надо отладить всё на идеальной модели входного сигнала, где всё точно.

С какой разрядностью на входе и выходе работает функция arm_cmplx_mag_q15()?

И ещё Вы пишете:
Если не будет выхода за предела разрядной сетки, сдвиг можно сделать хоть до, хоть после функции. Но Вы ещё не назвали разрядность, поэтому как будет известно - будем думать.

И сдвиг вправо, а не влево. Это ошибка в Ваших расчётах или просто опечатка? Если ошибка, то прошу все вышеприведённые цифры повторить заново после исправления.


Функция arm_cmplx_mag_q15() принимает на вход массив, состоящий 16 битных знаковых значений. В описании к функции написано следующее:

Цитата
The function implements 1.15 by 1.15 multiplications and finally output is converted into 2.14 format.


Я еще так и до конца и не разобрался, что значит форматы 1.15 и 2.14. Но из описания следует, что на вход подаем в формате 1.15, а на выходе имеем 2.14. Более того, входные значения у этой функции в комплексной форме, как я уже и писал - input[n] - действительная часть, input[n+1] - комплексная (n=0,2,4,6,8...). В эту функцию я отправляю массив только после БПФ, предварительно сдвинув каждый член на 10 бит ВЛЕВО, ошибки тут нет:

Код
for (u32 i=0; i<2048; i++) input[i] = input[i]<<10; //Насколько я понимаю -  это сдвиг ВЛЕВО, или я чего то конкретно упустил???


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

1) Сгенерировал 1024 точки синуса с частотой 50 Гц и амплитудой 600. Сигнал двуполярный один бит знак, остальное в дополнительном коде. Вот по такой схеме (с матлабом связываться не стал, проще было прям внутри моей программы его генерировать):

Код
    u32 inputShift=0;
    for (u32 i = 0; i < FFT_LEN; i++)
    {
        float tmp;
        tmp = TEST_AMPL * sin( (TEST_FREQUENCY*2*PI/(FFT_LEN-1)) * i);
        int tmpInt;
        tmpInt = tmp;
        input[inputShift++] = tmpInt;
        input[inputShift++] = 0;
    }


Получилась вот такая картинка:

Нажмите для просмотра прикрепленного файла

2) Далее полученные массив отправил на обработку в DSP либу, получилась вот такая характеритика:

Нажмите для просмотра прикрепленного файла
Krys
Цитата(yanvasiij @ Aug 24 2015, 19:06) *
Функция arm_cmplx_mag_q15() массив, состоящий 16 битных знаковых значений. В описании к функции написано следующее:
Мы говорим про функцию, которая рассчитывает амплитуду частотных составляющих по их мнимой и действительной частям?


Цитата(yanvasiij @ Aug 24 2015, 19:06) *
Я еще так и до конца и не разобрался, что значит форматы 1.15 и 2.14. Но из описания следует, что на вход подаем в формате 1.15, а на выходе имеем 2.14.
Ну я Вам про форматы пытался выше объяснить. Если что - спрашивайте конкретные вопросы, поясню более детально.
В рамках конкретной функции такие форматы можно объяснить: если каждую квадратуру выхода БПФ считать пронормированным к 1, то амплитуда в корень из 2 будет больше 1, поэтому при сохранении входного формата возникнет переполнение, следовательно положение точки нужно сдвинуть.
Но на самом деле этого не требуется конкретно для БПФ, т.к. заранее известно, что комплексные числа - это вращающиеся вектора, которые не могут выходить за единичную окружность, тогда и модуль не может превышать 1.


Цитата(yanvasiij @ Aug 24 2015, 19:06) *
предварительно сдвинув каждый член на 10 бит ВЛЕВО, ошибки тут нет:
А надо-то вправо. БПФ усилил на гейн. Надо вернуть обратно.


Цитата(yanvasiij @ Aug 24 2015, 19:06) *
1) Сгенерировал 1024 точки синуса с частотой 50 Гц и амплитудой 600
А частота дискретизации тогда какая?


Цитата(yanvasiij @ Aug 24 2015, 19:06) *
Получилась вот такая картинка:
2) Далее полученные массив отправил на обработку в DSP либу, получилась вот такая характеритика:
По этим картинкам видно, что частота сигнала у Вас не кратна частоте дискретизации. Подайте кратную частоту, так гораздо удобнее сначала разбираться, почему БПФ не дал ожидаемую амплитуду.
yanvasiij
Цитата(Krys @ Aug 25 2015, 09:52) *
Мы говорим про функцию, которая рассчитывает амплитуду частотных составляющих по их мнимой и действительной частям?


Ну я Вам про форматы пытался выше объяснить. Если что - спрашивайте конкретные вопросы, поясню более детально.
В рамках конкретной функции такие форматы можно объяснить: если каждую квадратуру выхода БПФ считать пронормированным к 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
Цитата(yanvasiij @ Aug 25 2015, 13:42) *
По поводу форматов конкретный вопрос собственно такой: У меня есть знаковое целочисленное число - отсчеты моего АЦП. Нужно ли его преобразовать, чтобы получить формат 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 бит влево.


Цитата(yanvasiij @ Aug 25 2015, 13:42) *
По поводу гейна и сдвига. В том-то и дело, если сдвигать ВПРАВО, то в выходном массиве после вызова arm_cmplx_mag_q15() одни нули. Если не сдвигать, то картинка характеристики "обрезается", нет "нижних" маленьких "тычков" (гармоники с малой амплитудой).
Понятно, Вы методом научного тыка решили, что сдвигать всё же надо )))
Я придерживаюсь мнения, что сдвигать вообще не надо никуда, функция БПФ и так всё сделала за Вас, оно и логично по-идее с точки зрения разработки юзер-френдли интерфейса, чтобы пользователь не парил голову, сколько стадий бабочки прошло, ему нужен только готовый результат.
Однако вот тут мне непонятно:
Цитата(yanvasiij @ Aug 25 2015, 13:42) *
Если не сдвигать, то картинка характеристики "обрезается", нет "нижних" маленьких "тычков" (гармоники с малой амплитудой).
Откуда у Вас после сдвига влево спектр приобретает "нижние маленькие тычки", если без сдвига их нет? До сдвига эта информация должна быть уже заложена в несдвинутый сигнал. А если она заложена, то как она обрезается отсутствием каких-либо действий (случай без сдвига)? Не должно такого быть. Подумайте, может, где-то ошибка.


Цитата(yanvasiij @ Aug 25 2015, 13:42) *
Частота дискретизации 1024 Гц. Повторил эксперимент для частоты сгенерированного математически сигнала 32 Гц (т.е. 1024 кратно 32), амплитуда 600. Картинка исходного сигнала вот:
На выходе вот такая картинка:
Другое дело. Вот с этого сразу и надо было начинать.
К стати, по картинке спектра: левая и правая палки имеют разный уровень. Этого быть не должно. Уж не знаю, откуда взялось. Может, функция кривая, может, мнимая часть где-то пролезла (хотя почти не реально). Может, другие участники форума подскажут, откуда берётся этот эффект. Если не найдётся причина - предлагаю входные числа затолкать в матлаб один к одному, и над ними проделать fft(). Там должно быть одинаково.
yanvasiij
Цитата(Krys @ Aug 25 2015, 12:28) *
...

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

...


Спасибо большое, за объяснения форматов! Наконец-то стало понятно.

А по-поводу сдвига. Я же сдвигаю ПЕРЕД тем как ложить в функцию arm_cmplx_mag_q15(), а в ВЫХОДНОМ массиве уже смотрю есть тычки или нет. Если не сдвигать вообще ВХОДНЫЕ ДАННЫЕ, то и в выходном массиве ничего кроме нулей нет, если сдвигать недостаточно сильно, то в выходном массиве "срезается" нижняя часть характеристики.
Krys
Пожалуйста!

А во входном массиве (т.е. на выходе БПФ) есть что-то кроме нулей в тех точках, где нет палок? Можете показать вещественный спектр и мнимый спектр до сдвига и до выделения амплитуды?
Даже если Вы сдвинули вход, чтобы повысить разрешение, то после этого нужно сдвинуть обратно, чтобы сохранить правильный масштаб чисел. И тогда всё равно получите нули ))
prig
Цитата(Krys @ Aug 25 2015, 10:28) *
...
Теперь вернёмся к конкретно нашим цифрам: на входе 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
Цитата(prig @ Aug 25 2015, 16:00) *
Звиздец какое объяснение.
Осторожнее с эмоциями, здесь им не место )



Цитата(prig @ Aug 25 2015, 16:00) *
Я так понимаю, Вы пытаетесь донести вариант реализации умножения 1.15*1.15 -> 1.15 с использованием целочисленного умножителя.
Не обязательно умножения. К стати, с умножением всё сложнее. Я скорее пытался объяснить, что значит такой формат в виде двух чисел с точкой, и почему он всего лишь мысленный.


Цитата(prig @ Aug 25 2015, 16:00) *
Диапазон значений корня из суммы квадратов вещественной и мнимой части в формате 1.15 (диапазон каждой что-то около -0.99997...0.99997) будет выходить за границы диапазона -1...1. Т.е., в общем случае, для представления модуля комплексного числа в вышеуказанном формате 16-ти битным числом потребуется фомат 2.14 (примерно -1.99994...1.99994). С округлением, ессно. Как именно получается этот результат, вопрос десятый.
Я там выше подобный механизм расписывал для бабочки:
Цитата(Krys @ Aug 25 2015, 11:52) *
В рамках конкретной функции такие форматы можно объяснить: если каждую квадратуру выхода БПФ считать пронормированным к 1, то амплитуда в корень из 2 будет больше 1, поэтому при сохранении входного формата возникнет переполнение, следовательно положение точки нужно сдвинуть.
Но на самом деле этого не требуется конкретно для БПФ, т.к. заранее известно, что комплексные числа - это вращающиеся вектора, которые не могут выходить за единичную окружность, тогда и модуль не может превышать 1.



Цитата(prig @ Aug 25 2015, 16:00) *
С учётом того, что ТС собирается использовать нормировочный коэффициент (см. вопрос), использование формата 2.14 как 1.15 это всего лишь вопрос значения этого нормировочного коэффициента. Результат этой функции никуда сдвигать не надо.
А такое объяснение не понял теперь я )))


prig
Цитата(Krys @ Aug 25 2015, 12:47) *
Осторожнее с эмоциями, здесь им не место )



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


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



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


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

С умножением всё достаточно просто.
Достаточно осознать, что из себя представляют диапазоны соответствующих форматов с фикс.точкой и диапазоны результатов операций.
Вот несколько наиболее показательных примеров:
Порядок диапазона для беззнаковых аргументов 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 - это фактически одно и то же. Различать их позволяет как раз нормировочный коэффициент.
Krys
Цитата(prig @ Aug 25 2015, 18:22) *
Осторожнее с изобретением объяснений типа "на пальцах".
Понятно: поперепирались по типу "ты дурак - сам дурак" )))
А как интересно Вы хотели, чтобы объяснялось? Просто 1:1 по учебнику? Дак ТС и сам его может прочитать. Ценность - в объяснении "своими словами".


Цитата(prig @ Aug 25 2015, 18:22) *
Осторожнее с изобретением объяснений типа "на пальцах". Можно попасть впросак.
А Вы этого сильно боитесь? Я не возражаю, если меня поправит более опытный коллега.


Цитата(prig @ Aug 25 2015, 18:22) *
С умножением всё достаточно просто.
Достаточно осознать, что из себя представляют диапазоны соответствующих форматов с фикс.точкой и диапазоны результатов операций.
Вот несколько наиболее показательных примеров:
Порядок диапазона для беззнаковых аргументов 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 и потом округлить.
Тут теперь уже Ваши объяснения очень тяжело понять, так что даже дискутировать не могу. Единственно результат сдвигать надо вправо в последнем предложении, а округлить надо до сдвига.
К стати: вопрос на засыпку для разминки мозгов: почему при умножении двух 16-разрядных знаковых чисел в допкоде результат получается 32-разрядным, когда по казалось бы очевидной логике должен быть 31 разряд? )))


Цитата(prig @ Aug 25 2015, 18:22) *
С БПФ у Вас как раз незадача и вышла.
Давайте, учите меня писать БПФ ))) Я реализовал БПФ для сверхдлинных последовательностей )))


Цитата(prig @ Aug 25 2015, 18:22) *
Если использовать арифметику 1.15*1.15 -> 1.15 и 1.15+-1.15 -> 1.15 , возможность переполнения при вычислении БПФ гарантирована.
Именно из-за этого используется 1 или 2 "защитных" бита для радикса-2 или радикса-4 соответстенно.
А я разве не про это же написал?
Цитата(Krys @ Aug 24 2015, 11:00) *
Бабочка выполняет умножение входного числа (отсчёта синуса) на поворачивающий множитель. Поворачивающий множитель по определению не превышает 1, поэтому его формат тоже выбран 1.15 (хотя, это чуть-чуть некорректно, и в своей реализации я выбрал честный формат). После умножения бабочка делает сложение. Отсюда разрядность должна увеличиться на 1 бит. Разрядная сетка у нас фиксирована 16 битами, поэтому приходится отбросить младший бит дробной части, чтобы не словить переполнение. Отсюда получается результирующий формат на выходе бабочки 2.14. Бабочка выполняется над массивом данных 1024 сэмпла 10 раз, поэтому в результате точка съезжает аж на 10 разрядов вправо, про это говорят "гейн БПФ".



Цитата(prig @ Aug 25 2015, 18:22) *
Т.е. весь входной массив нормируется до необходимых значений. Нормировочный коэффициент запоминается.
Зачем что-то запоминать, если заранее из теории БПФ известно, что бабочка делается 10 раз, и при каждой делается сдвиг на 1 разряд? Можно просто на выходе сделать обратный сдвиг на сразу известную константу.


Цитата(prig @ Aug 25 2015, 18:22) *
Кроме того, на таком формате обычно используется групповая плавающая точка. Т.е., постадийная нормировка по-необходимости. Нормировочный к-т обновляется, ессно.
С промежуточной арифметикой 2.30+-2.30 -> 2.30 всё немного по-другому, но постадийная нормировка нормировка всё равно обычно используется.

Нажмите для просмотра прикрепленного файла

На приведённом рисунке (моделировал когда-то для своей задачи) сиреневым (power) показана степень того самого масштабирующего коэффициента, про который Вы говорили, в зависимости от стадии обработки. Для блочной плавающей точки. Видно, что, начиная с некоторой стадии, степень возрастает на каждой стадии. А "задержка" начала возрастания степени объясняется тем, что входной сигнал был не на полную шкалу. Если бы на входе была проделана нормализация, то ступеньки бы пошли сразу. Таким образом, из рисунка видно, что блочная плавающая точка не даёт никакого выигрыша для БПФ по сравнению с фиксированной, безусловно устанавливаемой для каждой стадии, т.к. всё равно ступеньки идут на каждой стадии.
Krys
что-то ТС пропал... у Вас всё заработало?
prig
Цитата(Krys @ Aug 25 2015, 07:52) *
...
Ну я Вам про форматы пытался выше объяснить. Если что - спрашивайте конкретные вопросы, поясню более детально.
В рамках конкретной функции такие форматы можно объяснить: если каждую квадратуру выхода БПФ считать пронормированным к 1, то амплитуда в корень из 2 будет больше 1, поэтому при сохранении входного формата возникнет переполнение, следовательно положение точки нужно сдвинуть.
Но на самом деле этого не требуется конкретно для БПФ, т.к. заранее известно, что комплексные числа - это вращающиеся вектора, которые не могут выходить за единичную окружность, тогда и модуль не может превышать 1.
...


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

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


Цитата(Krys @ Aug 26 2015, 06:50) *
...
Тут теперь уже Ваши объяснения очень тяжело понять, так что даже дискутировать не могу. Единственно результат сдвигать надо вправо в последнем предложении, а округлить надо до сдвига.
К стати: вопрос на засыпку для разминки мозгов: почему при умножении двух 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), можно говорить о диапазоне значений результата или его порядке, а также о формате результата и диапазоне значений этого формата. И в этой части всё очевидно.


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


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

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



П.С. Судя по всему, Вы всё время пытаетесь интерпретировать единственный частный случай реализации целочисленного БПФ.
Это не есть правильно. Таки, лучше подучить матчасть для общего случая.
Krys
Цитата(prig @ Aug 27 2015, 17:59) *
Цитата(Krys @ Aug 25 2015, 11:52) *

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

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


Цитата(prig @ Aug 27 2015, 17:59) *
В зависимости от используемых целочисленных операций, сдвиг (нормализация) может осуществляться в процессе вычисления "бабочки" БПФ.
Приведите пожалуйста пример, где модуль выхода БПФ больше максимально возможного размаха его квадратур. Для векторов, вращающихся вокруг единичной окружности это невозможно.


Цитата(prig @ Aug 27 2015, 17:59) *
Вероятно, Вы пыталисьсь обобщить частный случай реализации целочисленного БПФ и сделали ошибку.
Я привёл условия типа "если каждую квадратуру выхода БПФ считать пронормированным к 1". Я согласен, что это для частного случая. Но для него я не сделал ошибку.


Цитата(prig @ Aug 27 2015, 17:59) *
Т.е., два старших бита результата умножения 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


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

Цитата(prig @ Aug 27 2015, 17:59) *
Аналогично, заданный Вами вопрос о разрядности тоже оказался не слишком корректен.
Ещё как корректен, просто это нужно прощупать "на своей шкуре".

Цитата(prig @ Aug 27 2015, 17:59) *
Для такой детерминированной функции, как целочисленное умножение со знаком (16.0*16.0 -> 32.0), можно говорить о диапазоне значений результата или его порядке, а также о формате результата и диапазоне значений этого формата. И в этой части всё очевидно.
Очевидно, но первое, что очевидно, оказывается неверным. Вопрос на засыпку был "почему при умножении двух 16-разрядных знаковых чисел в допкоде результат получается 32-разрядным, когда по казалось бы очевидной логике должен быть 31 разряд"? Очевидно это 31 разряд, а не 32. Потому что при умножении полезные данные сосредоточены в младших 15 битах операндов, следовательно 15 * 15 получаем 30 и 1 бит знака. 31 разряд в результате. Откуда 32?


Цитата(prig @ Aug 27 2015, 17:59) *
- Вероятно, Вы моделировали с табличным синусом или чем-то простеньким.
Согласен, тут я сделал свой вывод только для своего сигнала. Он был реальный, но с сильным преобладанием несущей, на фоне которой нужно было обнаружить отражения с допплером.


Цитата(prig @ Aug 27 2015, 17:59) *
П.С. Судя по всему, Вы всё время пытаетесь интерпретировать единственный частный случай реализации целочисленного БПФ.
Это не есть правильно. Таки, лучше подучить матчасть для общего случая.
А Вы уверены, что смогли охватить этот самый общий случай? Он бесконечен. Все случаи учесть невозможно. Так-то хорошо все утверждения парировать, что это не годится, в общем случае работать не обязано. Универсальный козырь.
Krys
Цитата(Krys @ Aug 25 2015, 11:52) *
А надо-то вправо. БПФ усилил на гейн. Надо вернуть обратно.
Я извиняюсь, запутался. Действительно, здесь надо влево. БПФ усилил, но, чтобы влезть в разрядность сдвинул вправо с потерей младших битов. А чтобы получить число в его реальном масштабе, нужно обратно сдвинуть влево.
prig
Цитата(Krys @ Aug 28 2015, 09:07) *
Я извиняюсь, запутался.
...


Ну, на самом деле, в целочисленной арифметике много путаницы вообще, и в терминологии в частности. Плюс, масса мелочей, которые можно не заметить. Так что немудрено. Со мной тоже случалсь, иногда с последствиями. К собственно предмету FFT это имеет только косвенное отношение, но жизнь подсластить может.

Мне в своё время пришлось довольно долго объяснять заказчику, почему для их задачи лучше использовать радикс-2, а не 4, а ещё лучше соскочить с 16-бит АДи на 32-бит ТИ и избавиться от нескольких неудачных решений. Уговорил. Но с терминологией и мелочами пришлось разбираться серьёзно и, опять же, доносить до заказчика. При том, что заказчик с ЦОС не по-наслышке знаком, и вообще весьма крут и хорошо известен в узких кругах.

Ну, да ладно. Базар мы слегка не по-делу устроили, но кому-то он может оказаться полезным.
Объясняльщики мы оба никакие, но суть возможных неприятностей с целочисленными БПФ наверное донесли. И то ладно.


А вот если по сабжу, то возникает банальный вопрос. Зачем ТС морочится с 16-бит БПФ на F4?
Ну, разве что память экономит, но не факт, что это правильно в перспективе дальнейшего развития.
Krys
куда делся ТС? Мы уже тут со скуки похоливарить успели, а его всё нет )))
yanvasiij
Цитата(Krys @ Aug 31 2015, 13:21) *
куда делся ТС? Мы уже тут со скуки похоливарить успели, а его всё нет )))


Я сильно извиняюсь, за то, что пропал. Меня отправили в командировку, не могу сейчас ответить по существу. Как только вернусь покажу, что у меня получилось.
alex2103
Цитата(yanvasiij @ Sep 1 2015, 11:44) *
Я сильно извиняюсь, за то, что пропал. Меня отправили в командировку, не могу сейчас ответить по существу. Как только вернусь покажу, что у меня получилось.

Расскажите чем все закончилось. Куда, что и когда надо сдвигать чтоб в итоге получить правильные значения амплитуды гармоник?
yanvasiij
Цитата(alex2103 @ Feb 9 2016, 17:35) *
Расскажите чем все закончилось. Куда, что и когда надо сдвигать чтоб в итоге получить правильные значения амплитуды гармоник?


У меня все красиво заработало. Спасибо всем кто помогал и принимал участие! Ниже привожу пример на базе все той же dsp библиотеки. Делаю следующее сначала генерирую сигнал, потом загоняю его в FFT:

CODE

#define FFT_LEN 1024
#define TEST_FREQUENCY 31.969f
#define TEST_AMPL 600.0f

#define TEST_AMPL3 500.0f
#define TEST_FREQUENCY3 64

#define TEST_AMPL4 400.0f
#define TEST_FREQUENCY4 72

#define TEST_AMPL5 300.0f
#define TEST_FREQUENCY5 96

#define TEST_AMPL6 200.0f
#define TEST_FREQUENCY6 128

#define TEST_AMPL2 100.0f
#define TEST_FREQUENCY2 160

...
for (u32 i = 0; i < FFT_LEN; i++)
{
float tmp;
tmp = TEST_AMPL * sin( (TEST_FREQUENCY*2*PI/(FFT_LEN-1)) * i);
tmp += TEST_AMPL2 * sin( (TEST_FREQUENCY2*2*PI/(FFT_LEN-1)) * i);
tmp += TEST_AMPL3 * sin( (TEST_FREQUENCY3*2*PI/(FFT_LEN-1)) * i);
tmp += TEST_AMPL4 * sin( (TEST_FREQUENCY4*2*PI/(FFT_LEN-1)) * i);
tmp += TEST_AMPL5 * sin( (TEST_FREQUENCY5*2*PI/(FFT_LEN-1)) * i);
tmp += TEST_AMPL6 * sin( (TEST_FREQUENCY6*2*PI/(FFT_LEN-1)) * i);
tmp += 600;
int tmpInt;
tmpInt = tmp;
fInput[inputShift] = tmp;
fInput[inputShift+1] = 0;
input[inputShift++] = tmpInt;
input[inputShift++] = 0;
}
for (u32 i = 0; i < FFT_LEN*2; i++)
input[i] = fInput[i];

arm_cfft_f32(&arm_cfft_sR_f32_len1024, fInput, 0, 1);
arm_cmplx_mag_f32 (fInput, fOutput, FFT_LEN);
fOutput[0] = (fOutput[0]*NORMIRATION_COEF)/2;
output[0] = fOutput[0];
for (u32 k=1; k<FFT_LEN; k++)
{
fOutput[k] = fOutput[k]*NORMIRATION_COEF;
output[k] = fOutput[k];
}
...


Если нужно я могу конечно выложить весь проект, но там сделано на базе STMCube, так что очень много всего лишнего и весит при этом такой проект невменяемо много.
Apparatchik
Цитата(yanvasiij @ Feb 15 2016, 06:22) *
Если нужно я могу конечно выложить весь проект, но там сделано на базе STMCube, так что очень много всего лишнего и весит при этом такой проект невменяемо много.

Очень хотелось бы увидеть более полную картину, где уже используется АЦП. Если не затруднит выложите пожалуйста.
yanvasiij
Выложу в понедельник - пример у меня на работе
Apparatchik
Буду ждать
yanvasiij
Приношу извинения за задержку, не было времени нормально все проверить. Выкладываю как есть, там все вперемешку, это даже не рабочий проект, а так просто попробовали.
Нажмите для просмотра прикрепленного файла
Для просмотра полной версии этой страницы, пожалуйста, пройдите по ссылке.
Invision Power Board © 2001-2024 Invision Power Services, Inc.