реклама на сайте
подробности

 
 
4 страниц V   1 2 3 > »   
Reply to this topicStart new topic
> FFT stm32f407 лишние гармоники в выходном массиве
yanvasiij
сообщение Aug 19 2015, 10:48
Сообщение #1


Местный
***

Группа: Свой
Сообщений: 315
Регистрация: 23-12-11
Из: Уфа
Пользователь №: 69 041



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

Процессор 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[])

Прикрепленное изображение
Go to the top of the page
 
+Quote Post
Krys
сообщение Aug 21 2015, 06:23
Сообщение #2


Профессионал
*****

Группа: Свой
Сообщений: 1 994
Регистрация: 17-01-06
Из: Томск, Россия
Пользователь №: 13 271



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


--------------------
Зная себе цену, нужно ещё и пользоваться спросом...
Go to the top of the page
 
+Quote Post
prig
сообщение Aug 21 2015, 07:26
Сообщение #3


Знающий
****

Группа: Свой
Сообщений: 717
Регистрация: 30-01-08
Из: СПб
Пользователь №: 34 595



Цитата(yanvasiij @ Aug 19 2015, 13:48) *
...
Однако в результате output[0] всегда равен нулю, в output[1] некоторое значение ...
...
Вдобавок нет нет да и появятся непонятные значения в разных гармониках (то output[200], то output[511] и т.п. совершенно случайным образом).
...

- Для "фирменных" библиотек опыт показывает, что если правильные данные данные в правильном формате правильно размещены во входном массиве, результат соответствует теории с точностью до эффектов насыщения и т.д. Намёк ясен?
- О спектре оконной функции у Вас есть представление? Если нет, то Вы явно не с того начали.
Go to the top of the page
 
+Quote Post
Krys
сообщение Aug 21 2015, 09:53
Сообщение #4


Профессионал
*****

Группа: Свой
Сообщений: 1 994
Регистрация: 17-01-06
Из: Томск, Россия
Пользователь №: 13 271



Не умничайте )) Напишите, что конкретно


--------------------
Зная себе цену, нужно ещё и пользоваться спросом...
Go to the top of the page
 
+Quote Post
Alex11
сообщение Aug 21 2015, 10:01
Сообщение #5


Профессионал
*****

Группа: Свой
Сообщений: 1 986
Регистрация: 23-10-04
Из: С-Петербург
Пользователь №: 965



Конкретно - если не наложить окно на измеренный сигнал перед Фурье, то можно узнать о себе много нового и интересного в зависимости от погоды в Лондоне и др. внешних факторов.
Go to the top of the page
 
+Quote Post
Lmx2315
сообщение Aug 21 2015, 10:54
Сообщение #6


отэц
*****

Группа: Свой
Сообщений: 1 472
Регистрация: 18-09-05
Из: Москва
Пользователь №: 8 684



Цитата(Alex11 @ Aug 21 2015, 14:01) *
Конкретно - если не наложить окно на измеренный сигнал перед Фурье, то можно узнать о себе много нового и интересного в зависимости от погоды в Лондоне и др. внешних факторов.

..по умолчанию окно всегда наложено - прямоугольное, что как бы ничего не объясняет.
з.ы.
я бы сделал 1024 отчёта табличного синуса и посмотрел как он выглядит, только синус и всё, потом он и постоянка (искусственная const).
Потом два синуса, потом синус на рабочей частоте.


--------------------
"..не нравятся мои выборы? ..приходите в мой суд."
Узурпатор П.
Go to the top of the page
 
+Quote Post
Krys
сообщение Aug 21 2015, 11:06
Сообщение #7


Профессионал
*****

Группа: Свой
Сообщений: 1 994
Регистрация: 17-01-06
Из: Томск, Россия
Пользователь №: 13 271



Ну я что и предлагал - тесты на заранее подготовленных данных, а не на "живых" ))


--------------------
Зная себе цену, нужно ещё и пользоваться спросом...
Go to the top of the page
 
+Quote Post
yanvasiij
сообщение Aug 21 2015, 11:43
Сообщение #8


Местный
***

Группа: Свой
Сообщений: 315
Регистрация: 23-12-11
Из: Уфа
Пользователь №: 69 041



Цитата(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) *
- Для "фирменных" библиотек опыт показывает, что если правильные данные данные в правильном формате правильно размещены во входном массиве, результат соответствует теории с точностью до эффектов насыщения и т.д. Намёк ясен?
- О спектре оконной функции у Вас есть представление? Если нет, то Вы явно не с того начали.


- А я и не выставлял претензий к "фирменным" библиотекам и в их работоспособности не сомневаюсь. Вопрос вообщем в том и состоит, правильно ли я размещаю данные и в том ли формате? Потому и привел код и попытался подробно объяснить, что я делаю, может я где ошибся или чего неправильно понял. Вопрос именно в этом и был.
- Если речь идет спектре сигнала ограниченного во времени, и вообще о том как должен выглядеть спектр на выходе с ДПФ, то да у меня есть об этом представление.
Go to the top of the page
 
+Quote Post
Alex11
сообщение Aug 21 2015, 11:44
Сообщение #9


Профессионал
*****

Группа: Свой
Сообщений: 1 986
Регистрация: 23-10-04
Из: С-Петербург
Пользователь №: 965



Может быть мне следовало изъясняться тщательнЕе. Естественно, прямоугольное окно есть всегда, только оно не обеспечивает отсутствие разрывов функции и ее производной на краях окна. В результате, при периодичности сигнала не кратной размеру окна, получаем кучу артефактов. На живом сигнале это бывает всегда в той или иной степени.
Go to the top of the page
 
+Quote Post
yanvasiij
сообщение Aug 21 2015, 11:52
Сообщение #10


Местный
***

Группа: Свой
Сообщений: 315
Регистрация: 23-12-11
Из: Уфа
Пользователь №: 69 041



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


Я, если честно не нашел ни одного упоминая о том, использует ли DSP stm32f4xx библиотека оконное сглаживание. Выходит, что если нет, то подобные артефакты могут появится? Если да, то как это можно проверить? Сгенерировать сигнал не кратный размеру окна?
Go to the top of the page
 
+Quote Post
Lmx2315
сообщение Aug 21 2015, 11:57
Сообщение #11


отэц
*****

Группа: Свой
Сообщений: 1 472
Регистрация: 18-09-05
Из: Москва
Пользователь №: 8 684



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

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


--------------------
"..не нравятся мои выборы? ..приходите в мой суд."
Узурпатор П.
Go to the top of the page
 
+Quote Post
yanvasiij
сообщение Aug 21 2015, 12:04
Сообщение #12


Местный
***

Группа: Свой
Сообщений: 315
Регистрация: 23-12-11
Из: Уфа
Пользователь №: 69 041



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


Да, вот только значение в output[50] (тот пик, что на картинке) равно 10543, а у меня АЦП 12 бит, т.е. максимум 4095. Как так?
Go to the top of the page
 
+Quote Post
prig
сообщение Aug 21 2015, 12:23
Сообщение #13


Знающий
****

Группа: Свой
Сообщений: 717
Регистрация: 30-01-08
Из: СПб
Пользователь №: 34 595



Цитата(Lmx2315 @ Aug 21 2015, 13:54) *
..по умолчанию окно всегда наложено - прямоугольное, что как бы ничего не объясняет.
з.ы.
я бы сделал 1024 отчёта табличного синуса и посмотрел как он выглядит, только синус и всё, потом он и постоянка (искусственная const).
Потом два синуса, потом синус на рабочей частоте.


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

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

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

И не найдёте. Оконная функция - это очень простая операция (или отсутствие оной в случае прямоугольного окна).
Выбор оконной функции далеко не элементарен и полностью лежит на совести разработчика.
Go to the top of the page
 
+Quote Post
Lmx2315
сообщение Aug 21 2015, 12:49
Сообщение #14


отэц
*****

Группа: Свой
Сообщений: 1 472
Регистрация: 18-09-05
Из: Москва
Пользователь №: 8 684



Цитата(yanvasiij @ Aug 21 2015, 16:04) *
Да, вот только значение в output[50] (тот пик, что на картинке) равно 10543, а у меня АЦП 12 бит, т.е. максимум 4095. Как так?

..ну и что? интерес практический или академический?
ваши 4095 были во временной области, а вы уже смотрите частотную , где совсем другие размерности.
Если вам надо сравнивать сигналы между собой - отнормируйте спектр по известному сигналу.
Подайте виртуальный синус размером с дин. диапазон вашего АЦП и посмотрите сколько попугаев он будет занимать в частотной области.
Кстати, тут уже будут сильно сказываться применяемые вами окна.


--------------------
"..не нравятся мои выборы? ..приходите в мой суд."
Узурпатор П.
Go to the top of the page
 
+Quote Post
Alex11
сообщение Aug 21 2015, 16:03
Сообщение #15


Профессионал
*****

Группа: Свой
Сообщений: 1 986
Регистрация: 23-10-04
Из: С-Петербург
Пользователь №: 965



Кроме всего изложенного выше, в спектре для определения амплитуды исходного сигнала нельзя использовать одиночный отсчет. Нужно брать квадратный корень из суммы квадратов отсчетов, попадающих в полосу сигнала. Ширина этой полосы определяется оконной функцией и требуемой точностью.
Go to the top of the page
 
+Quote Post

4 страниц V   1 2 3 > » 
Reply to this topicStart new topic
1 чел. читают эту тему (гостей: 1, скрытых пользователей: 0)
Пользователей: 0

 


RSS Текстовая версия Сейчас: 28th May 2017 - 03:01
Рейтинг@Mail.ru


Страница сгенерированна за 0.01688 секунд с 7
ELECTRONIX ©2004-2016