БПФ на 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

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