БПФ на LESO4 с использованием libFFTW

Рязанов Илья А., ryasanov@gmail.com

Когда-то давно основной задачей любого осциллографа было фиксация формы сигнала. И действительно, базовый функционал определил название, ведь дословно, осциллограф означает пишущий колебания (осцилло – качаюсь, граф – пишу ). Но современный осциллограф давно уже больше, чем "рисователь" сигнала, теперь это комплексный измерительный прибор, который помимо формы сигнала измеряет массу других связанных параметров: амплитудное, среднее, среднеквадратическое, средневыпрямленное напряжение, частота сигнала, всевозможные коэффициенты, плотность распредения и многое другое. Но особым пунктом всегда стоит спектр сигнала. Думаю, никому не нужно объяснять, какую роль играет спектральный анализ в современной метрологии. Потому при разработке анализатора сигналов LESO4 мы особое внимание уделили возможности исследовать спектр сигнала. Однако аппаратно в анализаторе сигналов LESO4, не реализована возможность расчета БПФ. Прибор по сути своей представляет собой умный многоканальный АЦП с аналоговой обвязкой. Используя API функции, можно сконфигурировать аналоговый измерительный тракт (открытый/закрытый вход, входной делитель/усилитель) и получить выборку сигнала при требуемой эквивалентной частоте дискретизации. Расчет спектра целиком и полностью ложится на программное обеспечение управляющей программы.

Получение спектра дискретного сигнала осуществляется дискретным преобразованием Фурье (ДПФ). Недостаток ДПФ состоит в очень низком быстродействии, поэтому на практике используют разновидность ДПФ - быстрое преобразование Фурье (БПФ), которое в зарубежной литературе называется - Fast Fourier Transform (FFT). В целом, алгоритм БПФ не является секретом и доступен в любом учебнике по цифровой обработки сигналов. При должном навыке программирования алгоритм можно написать самому, однако на изобретение велосипеда уйдет достаточно много времени, да и разобраться с алгоритмом придется на все 110%, потому как подводных камней предостаточно. В своей разработке мы пошли другим путем, решили использовать готовую библиотеку с уже написанным и полностью отлаженным алгоритмом. В этой статье попробую привести пример простенькой программы с использованием свободной библиотеки libFFTW (сайт проекта: http://www.fftw.org/), являющейся стандартом де-факто для многих программ в мире Open Source.

Особенности libFFTW:

Распространяется по лицензии GPL (да, на свое ПО мы тоже даем лицензию GPL).
Поддержка одномерного и многомерных БПФ преобразований.
Кроссплатформенность (можно скомпилировать на любой платформе, где есть компилятор С).
Поддержка многопоточности (POSIX thread и MPI).
Написана на C, теоретически можно использовать с любым языком программирования. Из коробки есть поддержка Fortran (Месье знает толк в извращениях).
Библиотека работает с комплексными и реальными числами.

Используя библиотеку, достаточно просто получить спектр сигнала, но визуализация спектра – самостоятельная и не простая задача, потому в нашем простом примере спектр используем для нахождения базовой основной частоты сигнала и в программе ограничимся текстовым режимом работы. Итак, задача: используя осциллограф LESO4 и библиотеку libFFTW, найти частоту сигнала.

Алгоритм будет простым. Сначала считаем с прибора выборку сигнала, пропустим через прямое БПФ. Затем в массиве полученных спектральных отсчетов находится элемент с наибольшим значением и запоминается его индекс. Полученный индекс умножим на шаг БПФ. Полученное значение и есть частота исследуемого сигнала. Это - основная частота, потому как в спектре может быть несколько гармоник. Нетрудно догадаться, что точность расчета частоты сигнала зависит от частоты дискретизации и количества выборки отсчетов. Для более сложных манипуляций со спектром не стоит забывать про окна сглаживания, которые придется реализовывать, но это совсем другая история.

Шаг БПФ определяется:

STEP_FFT = f_d / N,
где:
f_d - частота дискретизации сигнала;
N - Количество отсчетов;

Эти величины устанавливаются в приборе посредством вызова api функций:

leso4SetSamplFreq(sampe_frequency_25MHz);
leso4SetSamplesNum(8192);

В приборе 8192 – максимальный размер выборки. Тогда шаг дискретизации будет равен почти 3 кГц:

STEP_FFT = f_d / N = 25000000 / 8192 = 3051,7578125 Гц

Описание используемых функций библиотеки libFFTW

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

fftw_plan fftw_plan_dft_1d(int n, fftw_complex *in, fftw_complex *out, int sign, unsigned flags);

Создание плана одномерного БПФ преобразования. Функция принимает следующие параметры:

n – сколько отсчетов будем преобразовывать, размер выборки;
in – указатель на входной массив отсчетов;
out – указатель на массив выходных отсчетов;
sign – направления преобразования (FFTW_FORWARD — прямое БПФ преобразование, FFTW_BACKWARD – обратное БПФ преобразование, flags – дополнительные флаги преобразования);

Итак, "все идет по плану". Запускаем преобразования в соответствии с планом plan:

void fftw_execute(const fftw_plan plan);

Мы же хорошие программисты, прибираемся за собой: освобождаем память после преобразования. Функции передается дескриптор преобразования plan:

void fftw_destroy_plan(fftw_plan plan);

Тип данных fftw_complex, представляет собой массив из двух элементов

typedef double fftw_complex[2];

Первый элемент - реальная часть, второй - мнимая.

Ниже приведен листинг программы, рассчитывающий частоту сигнала:

Исходный код примера использования FFTW

// leso4_fft.c 
#include <stdio.h>
#include <leso4.h>
#include <fftw3.h>
#include <math.h>
 
#define DEVICE_NAME "LESO4.1_ER" // Дескриптор устройства
#define SAMPLES_NUM 8192 // Количество получаемых отсчетов
#define SAMPLES_FREQ 25000000 // Частота дискретизации
 
// Заполнить буфер sampleBufferPtr отсчетами, которые получены с устройства
// Так как нас не интересует фаза сигнала, то мнимые части делаем равными нулю
void getSample(fftw_complex *sampleBufferPtr, int sampleNum)
{
	int i = 0;
	double *complexDigit;
	double *deviceSample = (double*)leso4GetData(NULL, CHANNEL_A);
 
	for(i = 0; i < sampleNum; i++)
	{
		complexDigit = (double*)&sampleBufferPtr[i];
		complexDigit[0] = deviceSample[i];	// re - реальная часть
		complexDigit[1] = 0;			// im - мнимая часть
	}
}
 
// Расчет частоты по спектральным отсчетам
double calculateFreq(fftw_complex *sampleBufferPtr, int sampleNum)
{
	int i = 0;
	double *complexDigit;
	int maxIndex = 0;
	double maxDigit = 0;
 
	// Поиск максимального значения и запись его индекса в переменную maxIndex, тупым перебором
	for(i = 0; i < sampleNum; i++)
	{
		complexDigit = (double*)&sampleBufferPtr[i];
                // Расчет модуля комплексного числа
		double realDigit = 4*sqrt(pow(complexDigit[0], 2) + pow(complexDigit[1], 2))/sampleNum; 
		if(maxDigit < realDigit)
		{
			maxDigit = realDigit;
			maxIndex = i;
		}
	}
 
	return maxIndex*(SAMPLES_FREQ/SAMPLES_NUM); // Возвратить частоту сигнала
}
 
int main()
{
	fftw_plan forward_fft_plan;
	fftw_complex sampleBuffer[SAMPLES_NUM]; // Массив комплексных чисел
 
        // План прямого БПФ
        // Выходные спектральные отсчеты записываются в тот же массив, который передается в качестве входных отсчетов
        // Содержимое массива sampleBuffer, после выполнение БПФ перезаписывается
	forward_fft_plan = fftw_plan_dft_1d(SAMPLES_NUM, &sampleBuffer[0], &sampleBuffer[0], FFTW_FORWARD, FFTW_MEASURE);
 
	if(leso4Open(DEVICE_NAME) == LESO4_OK) // Открыть устройство
	{
		leso4EnableChannel(CHANNEL_A); // Включить канал A
		printf("Device %s open\n", DEVICE_NAME);
                printf("Version: %d\n", leso4GetVersion());
 
		leso4SetSamplFreq(sampe_frequency_25MHz); // Частота дискретизации АЦП 25 МГц
		leso4SetSamplesNum(SAMPLES_NUM); // Установить количество получаемых отсчетов равным 8192
		leso4SetInputMode(CHANNEL_A, input_mode_close); // Настроить на канале A закрытый вход
		leso4SetAmplitudeScan(CHANNEL_A, max_voltage_20V); // Установить аттенюатор канала А, на максимальное входное напряжение 20 В
 
		int rxBytes = leso4ReadFIFO(); // Прочитать данные с устройства
		printf("Read bytes: %d\n", rxBytes);
 
		getSample(&sampleBuffer[0], SAMPLES_NUM); // Заполнить массив отсчетами 
		fftw_execute(forward_fft_plan); // Выполнить план БПФ
 
                // Частота рассчитывается только для первой половины массива, потому-что во второй половине 
                // находятся те же отсчеты, что и в первой (см. свойства дискретного преобразования Фурье)
		double signalFreq = calculateFreq(&sampleBuffer[0], SAMPLES_NUM/2); 
 
		printf("Signal freq: %f Hz", signalFreq); // Вывести частоту сигнала на консоль
		leso4CloseDevice(); // Закрыть устройство
	}
	else // Девайс не получилось открыть
	{
		if(leso4GetLastErrorCode() == ERROR_CODE_FTDI_LIB_NOT_LOADED)
			printf("D2XX not install\n"); // Драйвер не установлен
		else
			printf("Error open device: %s\n", DEVICE_NAME); 
	}
	return 0;
}

Для компиляции программы должны быть установлены библиотеки libFFTW и libLESO4.
Компиляция в Linux или в Windows с помощью MinGW:

gcc -o leso leso4_fft.c -I{путь к заголовочным файлам} -L{путь к динамическим библиотекам} -lLESO4 -lfftw3-3

Как видно, если использовать готовые библиотеки, то ничего сложного в реализации ЦОС на ПК нет. Даже на языке программирования С. Этот прием расчета частоты я использовал в программе LESOscope для анализатора сигналов LESO4. Программная обработка отсчетов реального сигнала на компьютере очень удобна для разработки различных алгоритмов ЦОС в разных стендах и просто для изучения уже написанных алгоритмов.

Анализатор сигналов leso4, тот самый который используется в статье
Библиотека для работы с анализатором сигналов leso4
Библиотека FFTW для расчета БПФ
Кроссплатформенный виртуальный прибор, написанный на Qt, для анализатора сигналов leso4

16 января 2016

Добавить комментарий

Сообщение об ошибке

Новости ЛЭСО

Контактные лица

БПФ на LESO4 с использованием libFFTW

Описание используемых функций библиотеки libFFTW

Похожие статьи

Сообщение об ошибке

Вы здесь

Форма поиска

Новости ЛЭСО

Контактные лица

БПФ на LESO4 с использованием libFFTW

Описание используемых функций библиотеки libFFTW

Похожие статьи