Перейти к содержанию
    

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

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

 

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

 

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

 

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

 

post-69041-1439984938_thumb.png

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

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

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

...

Однако в результате output[0] всегда равен нулю, в output[1] некоторое значение ...

...

Вдобавок нет нет да и появятся непонятные значения в разных гармониках (то output[200], то output[511] и т.п. совершенно случайным образом).

...

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

- О спектре оконной функции у Вас есть представление? Если нет, то Вы явно не с того начали.

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Конкретно - если не наложить окно на измеренный сигнал перед Фурье, то можно узнать о себе много нового и интересного в зависимости от погоды в Лондоне и др. внешних факторов.

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Конкретно - если не наложить окно на измеренный сигнал перед Фурье, то можно узнать о себе много нового и интересного в зависимости от погоды в Лондоне и др. внешних факторов.

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

з.ы.

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

Потом два синуса, потом синус на рабочей частоте.

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

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

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Одной из причин может быть то, что подаваемая частота 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 бит влево. Сместил, совершенно не понимая зачем. В итоге код принял вот такой вид:

 

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

 


#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);
}

 

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

 

post-69041-1440156815_thumb.png

 

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

 

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

- О спектре оконной функции у Вас есть представление? Если нет, то Вы явно не с того начали.

 

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

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

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

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

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

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

 

Я, если честно не нашел ни одного упоминая о том, использует ли DSP stm32f4xx библиотека оконное сглаживание. Выходит, что если нет, то подобные артефакты могут появится? Если да, то как это можно проверить? Сгенерировать сигнал не кратный размеру окна?

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

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

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

думаю в выходном массиве точно не Дб, а линейные величины.

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

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

думаю в выходном массиве точно не Дб, а линейные величины.

 

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

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

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

з.ы.

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

Потом два синуса, потом синус на рабочей частоте.

 

- Если накладыватель не в курсе, считайте, что не наложено (последствия описаны Alex11).

- Если у ТС "выбивает" постоянку, что ему дадут табличные синусы?

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

М.б. битреверс там отдельной функцией, и он не приложился им в нужном месте.

 

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

 

Я, если честно не нашел ни одного упоминая о том, использует ли DSP stm32f4xx библиотека оконное сглаживание. Выходит, что если нет, то подобные артефакты могут появится? Если да, то как это можно проверить? Сгенерировать сигнал не кратный размеру окна?

И не найдёте. Оконная функция - это очень простая операция (или отсутствие оной в случае прямоугольного окна).

Выбор оконной функции далеко не элементарен и полностью лежит на совести разработчика.

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

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

..ну и что? интерес практический или академический?

ваши 4095 были во временной области, а вы уже смотрите частотную , где совсем другие размерности.

Если вам надо сравнивать сигналы между собой - отнормируйте спектр по известному сигналу.

Подайте виртуальный синус размером с дин. диапазон вашего АЦП и посмотрите сколько попугаев он будет занимать в частотной области.

Кстати, тут уже будут сильно сказываться применяемые вами окна.

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

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

Поделиться сообщением


Ссылка на сообщение
Поделиться на другие сайты

Присоединяйтесь к обсуждению

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

Гость
Ответить в этой теме...

×   Вставлено с форматированием.   Вставить как обычный текст

  Разрешено использовать не более 75 эмодзи.

×   Ваша ссылка была автоматически встроена.   Отображать как обычную ссылку

×   Ваш предыдущий контент был восстановлен.   Очистить редактор

×   Вы не можете вставлять изображения напрямую. Загружайте или вставляйте изображения по ссылке.

×
×
  • Создать...