Наблюдаем спектр сигнала при помощи D/QtE5 и Arduino

В этой статье мы покажем как сделать простой «анализатор спектра» почти в реальном времени своими руками и даже попробуем посмотреть некий сигнал, который запросто можно поймать на один медный провод (или даже самодельную антенну), присоединенный к плате Arduino.

И да, нечто подобное мы уже делали, но в этот раз будет уже иной уровень, и кроме того, можно вместо Arduino использовать даже iCEStick!

Для начала экспериментальной работы нам потребуются:

  • Плата Arduino Uno R3 (или ее аналог)
  • Медный провод или другая подручная антенна
  • Arduino IDE
  • Установленный компилятор D с менеджером пакетов dub
  • Библиотека serialport
  • Библиотека QtE5
  • Содинение с интернетом
  • Прямые руки

Создаем пустой проект с помощью dub:

dub init qt-receiver

Далее, создадим графический интерфейс, который будет выполнять отображение сигнала и его спектра мощности, который мы получим с помощью процедуры быстрого преобразования Фурье (но об этом чуть позже). Также, создадим удобные боксы для выбора подходящего для получения данных COM-порта (т.е. порта на который будет подключена плата) и для выбора нужной оконной функции, которая будет выполнять сглаживание рассчитанного спектра. На форме также разместим удобный слайдер, с помощью которого можно будет изменить ширину окна для преобразования Фурье (т.е. количество точек в процедуре этого преобразования), помимо этого возможно будет изменить этот параметр прямо во время работы процедуры получения данных.

В проекте будут использованы всего две кнопки с говорящими названиями Start и Stop, которые будут соответственно запускать и останавливать процедуру получения данных от микроконтроллера на плате, а небольшая панель статуса в нижней части окна программы будет показывать текущее состояние считывания данных.

Для создания интерфейса мы воспользуемся библиотекой QtE5: скопируем файл qte5.d в папку /source/ основного проекта и приступим к созданию интерфейса и не забудем про то, что у нас будет режим «почти реального времени», что означает помещение в программу и невизуального элемента из QtE5, а именно, таймера QTimer. Помимо таймера потребуются компоненты для отображения формы сигнала и спектральной картинки, которых нет в QtE5, что приводит к тому, что их нам придется написать самим, а пока выполним размещение элементов так, как будто эти компоненты у нас уже есть, для чего обозначим эти компоненты как QSignalBox (форма сигнала) и QSonogramBox (спектральная картина) и выполним размещение элементов в форме:

alias WindowType = QtE.WindowType;
alias normalWindow = WindowType.Window;

class MainForm : QWidget
{
	private
	{
		QHBoxLayout horizontalBox0, horizontalBox1, horizontalBox2;
		QVBoxLayout mainSizer, verticalBox0, verticalBox1, verticalBox2, verticalBox3,
					verticalBox4, verticalBox5, verticalBox7;
		QGroupBox group0, group1, group2, group3, group4, group5, group6, group7;
		QComboBox combo0, combo1;
		QSlider slider0;
		QLabel label0;
		QPushButton startButton, stopButton, overlayButton;
		QSignalBox area0;
		QSonogramBox area1;
		QStatusBar status0;
		QAction action0, action1, action2, action3;
		QTimer timer0;
	}

	private
	{
		USBDevice device;

		float[] data;
		WindowFunction windowFunction;
		uint windowSize;
	}

	this(QWidget parent, WindowType windowType)
	{
		super(parent, windowType);
		showMaximized;
		setWindowTitle("Digital Receiver Experiment");

		timer0 = new QTimer(this);
		timer0.setInterval(1);

		action0 = new QAction(this, &onTimerTick, aThis);
		connects(timer0, "timeout()", action0, "Slot()");

		mainSizer = new QVBoxLayout(null);

		horizontalBox0 = new QHBoxLayout(null);

		verticalBox0 = new QVBoxLayout(null);
		area0 = new QSignalBox(this);
		area0.saveThis(&area0);
		area0.setFixedHeight(80);
		area0.saveThis(&area0);

		verticalBox0.addWidget(area0);

		group0 = new QGroupBox(null);
		group0.setText("Signal data:");
		group0.setLayout(verticalBox0);
		group0.setFixedHeight(150);

		verticalBox1 = new QVBoxLayout(null);
		area1 = new QSonogramBox(this);
		area1.saveThis(&area1);

		verticalBox1.addWidget(area1);

		group1 = new QGroupBox(null);
		group1.setText("Sonogram:");
		group1.setLayout(verticalBox1);

		verticalBox3 = new QVBoxLayout(null);
		verticalBox3
			   .addWidget(group0)
			   .addWidget(group1);

		verticalBox4 = new QVBoxLayout(null);

		verticalBox5 = new QVBoxLayout(null);

		combo0 = new QComboBox(null);

		SerialPort
				 .listAvailable
				 .filter!(a => a.startsWith("/dev/ttyUSB"))
				 .enumerate(0)
				 .each!(a => combo0.addItem(a[1], a[0]));

		verticalBox5.addWidget(combo0);

		group3 = new QGroupBox(null);
		group3.setLayout(verticalBox5);
		group3.setText("COM port");
		group3.setMaximumHeight(70);

		horizontalBox2 = new QHBoxLayout(null);

		slider0 = new QSlider(null);
		slider0.setMinimum(2);
		slider0.setMaximum(9);
		slider0.setSingleStep(2);
		slider0.setSliderPosition(8);

		label0 = new QLabel(null);
		windowSize = 2 ^^ (slider0.value);
		label0.setText("<h3>%d</h3>".format(windowSize));

		horizontalBox2
				.addWidget(slider0)
				.addWidget(label0);

		group4 = new QGroupBox(null);
		group4.setLayout(horizontalBox2);
		group4.setText("Window size");
		group4.setMaximumHeight(70);


		horizontalBox1 = new QHBoxLayout(null);

		startButton = new QPushButton("Start", this);
		stopButton = new QPushButton("Stop", this);
		stopButton.setEnabled(false);

		horizontalBox1
					 .addWidget(startButton)
					 .addWidget(stopButton);

		action1 = new QAction(null, &onStartButton, aThis);
		action2 = new QAction(null, &onStopButton, aThis);
		action3 = new QAction(null, &onValueChanged, aThis);

  	        connects(startButton, "clicked()", action1, "Slot()");
		connects(stopButton, "clicked()", action2, "Slot()");
		connects(startButton, "clicked()", timer0, "start()");
		connects(stopButton, "clicked()", timer0, "stop()");
		connects(slider0, "sliderReleased()", action3, "Slot()");

		combo1 = new QComboBox(null);
		[
			"Rectangular",
			"Blackman",
			"Blackman(e)",
			"Blackman-Harris",
			"Blackman-Nutall",
			"Flat top",
			"Hamming",
			"Nutall",
			"Sine",
			"Hann"
		]
				 .enumerate(0)
				 .each!(a => combo1.addItem(a[1], a[0]));

		verticalBox7 = new QVBoxLayout(null);
		verticalBox7.addWidget(combo1);

		group7 = new QGroupBox(null);
		group7.setLayout(verticalBox7);
		group7.setText("Window function");
		group7.setMaximumHeight(70);


		verticalBox4
					.addWidget(group3)
					.addWidget(group7)
					.addWidget(group4)
					.addWidget(new QWidget(null))
					.addLayout(horizontalBox1);

		group6 = new QGroupBox(null);
		group6.setLayout(verticalBox4);
		group6.setText("Parameters");
		group6.setMaximumWidth(250);

		horizontalBox0
			.addLayout(verticalBox3)
			.addWidget(group6);

		status0 = new QStatusBar(null);
		status0.setMaximumHeight(15);
		status0.showMessage("Ready");

		mainSizer
				.addLayout(horizontalBox0)
				.addWidget(status0);

		setLayout(mainSizer);
	}
}

В этом длинном фрагменте кода мы производим объявление всех тех элементов, которые мы будем использовать в графическом интерфейсе, а также ряд определений которые позволят активно манипулировать данными, которые поступят с платы: device — объект COM-порта, windowFunction — перечисление для выбора оконной функции (об этом чуть ниже) и windowSize — целочисленная переменная для хранения размеров «окна». Дальше производим наследование от базового класса формы, устанавливаем максимальный размер окна с помощью метода showMaximized, создаем элементы интерфейса с попутным их размещением на форме.

Также, с помощью функционального подхода и унифицированного синтаксиса вызова функции заполняем бокс для списка имеющихся в системе COM-портов, осуществляя на ходу предварительную фильтрацию (не является необходимой, но очень удобна в среде Linux, если собираете под Windows, то необходимо убрать функцию filter). Бокс с выбором порта не единственный бокс, который требует заполнения во время старта программы, т.к. на форме есть бокс выбора оконной функции. Он заполняется названиями самых известных оконных функций, хорошо зарекомендовавших себя в сфере цифровой обработки сигналов (используемые окна в порядке перечисления их в программе: прямоугольное, Блекмана, расширенное Блекмана, Блекмана-Харриса, Блекмана-Наталла, «плоская шляпа», Хемминга, Наталла, синусоидальное и Ханна).

Теперь после размещения элементов интерфейса и подключения к ним нужных событий с помощью QAction и connects, опишем «переходники» для кнопок, таймера и слайдера:

extern (C)
{
    void onStartButton(MainForm* mainFormPointer)
    {
        (*mainFormPointer).runStart;
    }

    void onStopButton(MainForm* mainFormPointer)
    {
        (*mainFormPointer).runStop;
    }

    void onTimerTick(MainForm* mainFormPointer)
    {
        (*mainFormPointer).runTimer;
    }

    void onValueChanged(MainForm* mainFormPointer)
    {
        (*mainFormPointer).runChange;
    }
}

Оставим на некоторое время в покое почти готовую форму окна программы, сохранив ее в файле mainform.d и перейдем к описанию самой критической части нашей программы — описанию получения данных по USB посредством виртуального COM-порта, создаваемого чипом-контроллером на плате.

Не удивляйтесь такому резкому прыжку от дизайна графического интерфейса к работе с устройством на базе микроконтроллера, поскольку от того, как будет произведено описание взаимодействия компьютера с платой зависит функциональная структура всей остальной программы. Дело в том, что формат данных с микроконтроллера и способ их обработки влияет  на информационное наполнение двух виджетов QSignalBox и QSonogramBox, а также на вызов по таймеру определенных процедур подготовки и визуализации полученной информации.

Для начала работы по приему данных с USB порта, к которому будет подключена наша плата с Arduino, необходимо уточнить тот момент, что пересылаемыми данными будут результаты работы аналогово-цифрового преобразователя (АЦП). Этот преобразователь установлен на самой плате и является 10-разрядным (т.е. результатом преобразования аналогового сигнала является число в интервале от 0 до 1023 включительно) и именно этот момент несколько напрягает: встроенная функция Serial.write в Arduino позволяет передать только один байт, в который весь результат преобразования явно не влезет, а округлять значение до байта значит потерять точность результата…

К счастью, в нашем блоге уже описывалась идея, как можно это обойти: достаточно просто воспользоваться системой счисления по основанию 255, т.е. передавать два байта один из которых содержит результат деления на 255, а второй — только остаток от такого деления.

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

Пример модуля  приемника данных через COM-порт в файле usb.d:

module receiver.usb;

private
{
    import std.conv;
    import std.stdio;
    import std.string;

    import serialport;
}

class USBDevice
{
    private SerialPortBlk device;
    enum PORT_SPEED = 230_400;

    this(string name)
    {
        device = new SerialPortBlk(name, PORT_SPEED);
        device.config = SPConfig(PORT_SPEED, DataBits.data8, Parity.none, StopBits.one);
    }

    float read()
    {
        float result = 0.0f;

        try
        {
            ubyte[2] buffer;
            device.read(buffer);

            auto tmp = buffer[0] * 255 + buffer[1];

                if (tmp <=  1023)
                {
                    result = ((buffer[0] * 255 + buffer[1]) * 5.0f) / 1024.0f;
                }
                else
                {
                    result = -1.0f;
                }
        }
        catch (Throwable)
        {

        }

        return result;
    }

    void close()
    {
        device.close;
    }
}

В данном случае, мы используем несколько иную библиотеку для работы с COM-портом, чем в упоминавшейся статье, и потому название устройства и способ его конфигурации другие.

В этом примере мы создаем объект некоторого устройства, задаем для него имя порта и скорость (скорость у нас будет фиксированная, и составляет 230400 бит/с или 230400 бод), также устанавливаем 8 бит данных, без проверки четности, 1 стоп-бит. Метод read класса USBDevice предлагает надежный способ считывания данных, предлагая сборку значения полученного из двух байтов и вводя проверку для попадания данных в диапазон от 0 до 1023 (т.е. если в ходе сборки мы получаем значение большее, чем 1023, то это означает то, что мы приняли данные с USB неверно, а именно, что-то пропустили) — и если значение не попало в искомый диапазон, то результатом приема будет специальная величина -1.0, которая послужит нам в дальнейшем специальным маркером неудачно полученного значения.

Стоит заметить, что собранное значение умножается на 5.0, а затем делиться на 1024.0, что переводит значение из шкалы АЦП в значение напряжения на входе АЦП в вольтах (5.0 — максимальное напряжение на входе АЦП, 1024.0 — количество уровней напряжения).

Метод close используется для безопасного закрытия порта, который будет использоваться не только для закрытия порта, но и для остановки считывания данных в «реальном времени».

Теперь, напишем код передачи данных для платы Arduino, который почти не отличается от того, что мы упоминали ранее в одной из статей:

void setup()
{
    Serial.begin(230400);
}

void loop() {
    float value = analogRead(A0);

    byte low =  ((int) value) / 255;
    byte high = ((int) value) % 255;

    Serial.write(low);
    Serial.write(high);
}

После окончания работы с платой и кодом для устройства возвращаемся к описанию графического интерфейса, а именно двух основных виджетов: виджет визуализации сигнала (амплитуда) и виджет отображения его спектра. Описания обоих виджетов состоят из трех основных частей, которые функционально выполняют одну и ту же работу — графическое отображение некоторых данных. Итак, тремя составными частями виджета являются следующие методы: конструктор класса, сеттер полученных данных (т.е. метод, которым в виджет подается очередная порция данных) и метод отрисовки виджета, который запускает всю процедуру отображения самого виджета и вызывает внешнюю процедуру прорисовки данных.

Таким образом, виджеты отличаются лишь незначительными деталями в получении данных и процедурами визуализации данных.

Виджет отображения сигнала QSignalBox:

module receiver.qsignalbox;

private
{
    import qte5;

    import receiver.observer;

    extern (C)
    {
        void onPaintSignal(QSignalBox* signalBox, void* eventPointer, void* painterPointer)
        {
            (*signalBox).runPaint(eventPointer, painterPointer);
        }
    }

    void drawSignal(QPainter painter, ref float[] data, int w, int h)
    {
        float scaleY = h / 5.0f;

        foreach (x, y; data)
        {
            float X = x * (w / cast(float) data.length);
            float Y = h - (y * scaleY);

            painter.drawLine(cast(int) X, h, cast(int) X, cast(int) Y);
        }
    }
}

class QSignalBox : QWidget
{
    private
    {
        float[] data;
        bool isDrawing;
    }

    this(QWidget parent)
    {
        super(this);
        setStyleSheet("background: Black");
        setPaintEvent(&onPaintSignal, aThis());
    }

    void setData(ref float[] data)
    {
        this.data = data;
        this.isDrawing = true;
    }

    void runPaint(void* eventPointer, void* painterPointer)
    {
        QPainter painter = new QPainter('+', painterPointer);

        QColor color = new QColor(this);
        color.setRgb(
            cast(int) (255 * 0.1f),
            cast(int) (255 * 0.1f),
            cast(int) (255 * 0.8f),
            70
        );

        QPen pen = new QPen;
        pen.setColor(color);

        auto w = this.width;
        auto h = this.height;

        pen.setWidth(1);

        painter.setPen(pen);

        if (isDrawing)
        {
            isDrawing = false;

            drawSignal(painter, data, w, h);
            data = [];
        }

        painter.end;
    }
}

Конструктор класса, его структура, также как и «переходник» для него описывались в статье про создание своего виджета, и не представляет технического интереса. Метод setData устанавливает текущие данные во внутреннюю переменную и устанавливает в true внутреннюю переменную, разрешающую перерисовку. Если переменная isDrawing установлена, то происходит вызов процедуры drawSignal, в которую передаются в качестве параметров текущий QPainter, данные для отображения, длина и ширина виджета.

Процедура drawSignal вычисляет масштаб по оси Y (это масштаб в вольтах, а поскольку максимальное значение — 5 В, то происходит деление на эту величину, т.е. происходит масштабирование данных в диапазон от 0 до 5 В), далее происходит отображение каждого элемента данных через отображение его в масштабе по оси X (масштабируем сигнал по длине виджета, что определяется точно также как и для масштабирования оси Y, упомянутое выше) и по оси Y (поскольку, точка с координатами 0;0 — это верхний левый угол, то для отображения величины сигнала по оси Y, требуется передвинуть начало координат вниз и применить масштаб по оси Y).

Внутри процедуры runPaint выполняются уже самые обычные процедуры: создание QPainter, создание цвета, определение длины и ширины виджета, задание цвета для отображения сигнала и задание «карандаша» для рисования».

Для создания виджета отображения спектра в «реальном времени» (т.е. сонограммы), нам потребуется провести предварительную подготовку данных и определить структуры данных для хранения результатов обработки. Все дело в том, что мы будем изображать спектр сигнала с учетом того, что сигнал меняется во времени, поэтому необходимо создать такую структуру данных, которая бы хранила меняющийся во времени сигнал.

Смотрите, допустим мы получили данные по сигналу с USB и отобразили их, но этот же сигнал будет обрабатываться с помощью преобразования Фурье для получения его спектра мощности.

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

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

Сдвигающийся стек на базе обычного стека в файле structures.d:

module receiver.structures;

private
{
    import std.range;
}

// стек для хранения данных
class Stack(T)
{
    private:
        T[] elements;

    public:

        // добавление элемента в конец
        void push(T element)
        {
            elements ~= element;
        }

        // добавление ряда элементов
        void push(T[] elementsArray)
        {
            elements ~= elementsArray;
        }

        // удалить один элемент
        void pop()
        {
            --elements.length;
        }

        // количество элементов
        size_t length() const @property
        {
            return elements.length;
        }

        // вершина стека
        T top() @property
        {
            return elements[$-1];
        }

        // содержимое стека
        T[] content() @property
        {
            return elements;
        }

        // очистить стек
        void clear()
        {
            elements.length = 0;
        }

        // стек пустой ?
        bool empty() const @property
        {
            if (elements.length == 0)
            {
                return true;
            }
            else
            {
                return false;
            }
        }
}


class ShiftedStack(T, size_t M) : Stack!T
{
    private:
        T[M] emptyStack = T.init;

    public:
        this()
        {
            elements ~= emptyStack;
        }


        override
        {
            void push(T element)
            {
                elements ~= element;
                elements = elements[1..$];
            }


            void push(T[] elementsArray)
            {
                if (!elements.empty)
                {
                    elements ~= elementsArray;
                    elements = elements[0..M];
                }
            }

            void pop()
            {
                if (!elements.empty)
                {
                    super.pop();
                    elements ~= emptyStack;
                    elements = elements[0..M];
                }
            }

            size_t length() const @property
            {
                return M;
            }

            void clear()
            {
                elements = emptyStack;
            }
        }

        ref T get(size_t index)
        {
            return elements[index];
        }
}

Виджет QSonogramBox, построенный на базе сдвигающегося стека (файл qsonogrambox.d):

module receiver.qsonogrambox;

private
{
    import std.algorithm;
    import std.conv;
    import std.math : abs;
    import std.range;
    import std.string;
    import std.stdio;

    import receiver.structures;
    alias FLOW = float[];

    import qte5;

    extern (C)
    {
        void onPaintSonogram(QSonogramBox* sonogram, void* eventPointer, void* painterPointer)
        {
            (*sonogram).runPaint(eventPointer, painterPointer);
        }
    }

    void drawSonogramLine(QPainter painter, ref float[] data, int index)
    {
        auto MIN_N_MAX = data.fold!(min, max);
        enum COLOR_MIN = 0x02;
        enum COLOR_MAX = 0xff;
        float COLOR_STEP = (COLOR_MAX - COLOR_MIN) / (MIN_N_MAX[1] - MIN_N_MAX[0]);

        QColor color = new QColor(null);
        QPen pen = new QPen;
        pen.setWidth(2);

        foreach (x, y; data)
        {
            int COLOR_INTENSITY = cast(int) (COLOR_STEP * y);

            if (COLOR_INTENSITY < COLOR_MIN) { COLOR_INTENSITY = COLOR_MIN; } if (COLOR_INTENSITY > COLOR_MAX)
            {
                COLOR_INTENSITY = COLOR_MAX;
            }

            color.setRgb(
                cast(int) abs(0.8 * COLOR_INTENSITY),
                0,
                COLOR_MAX - COLOR_INTENSITY,
                COLOR_INTENSITY,
            );

            pen.setColor(color);

            painter.setPen(pen);

            painter.drawPoint(cast(int) x * 2, 512 - index);
            //painter.drawPoint((512 - index) * 2, cast(int) (x * 512.0f / data.length));
        }
    }
}

class QSonogramBox : QWidget
{
    private
    {
        ShiftedStack!(FLOW, 512) stack;
        bool isDrawing;
    }

    this(QWidget parent)
    {
        super(this);
        stack = new ShiftedStack!(FLOW, 512);
        setStyleSheet("background: Black");
        setPaintEvent(&onPaintSonogram, aThis);
    }

    void setData(ref float[] data)
    {
        stack.push(data);
        isDrawing = true;
    }

    void runPaint(void* eventPointer, void* painterPointer)
    {
        QPainter painter = new QPainter('+', painterPointer);

        if (isDrawing)
        {
            isDrawing = false;

            foreach (y; 0..512)
            {
                drawSonogramLine(painter, stack.get(y), y);
            }
        }

        painter.end;
   }
}

Как уже упоминалось, строение класса QSonogramBox похоже на строение QSignalBox, а именно, представлены конструктор класса, метод setData, метод runPaint. Внутри класса присутствует переменная stack, которая реализует сдвигающийся стек на 512 элементов. В эту переменную подставляются данные, которые будут получены из результатов преобразования Фурье (потерпите немного, чуть позже объясню, как), и именно этот стек используется для отрисовки сонограммы. В методе runPaint происходит отрисовка одной линии спектра, для чего используется функция drawSonogram, а сам runPaint с помощью foreach цикла проходит все линии, которые попали в стек.

Процедура drawSonogram используя индекс строки, текущий QPainter и данные, взятые из стека, рассчитывает минимальный и максимальный элемент в линии спектра, на основе полученных данных вычисляет цвет каждого конкретного элемента и осуществляет отображение каждого элемента данных в его цвете согласно полученному спектру мощности. Стоит упомянуть о том, что формула для вычисления цвета не фигурирует ни в одном из источников, и ее пришлось искать самостоятельно и не факт, что она действительно верна. Вся процедура отрисовки спектральных линий устроена достаточно элегантно, хотя и является одной из узких мест во всей программе.

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

В файле fft.d мы определяем сначала ряд окон для преобразований: сначала в виде перечисления WindowFunction, а затем в виде самостоятельных функций, каждая из которых принимает еще необработанные данные в виде изменяемого массива. После этого определяем метод applyWindow, принимающий набор данных и перечисление WindowFunction, который применяет выбранное окно (заданное окно является членом перечисления). Также мы определяем функцию, которая по строковому названию окна создаст член перечисления, что будет использовано нами для выбора окна из списка в боксе (данная функция есть практическое применение такого приема как фабричный метод, который обычно применяется для создания классов и является одним из паттерноа проектирования).

Также, определяем функцию получения спектра мощности makeFFT, которая сначала применяет окно к данным, делает преобразование Фурье (используя стандартную процедуру из std.numeric), вычисляет спектр мощности по соотвествующей формуле для кратковременного преобразования Фурье.

Внизу описания располагается функция, которая центрирует полученный спектр, располагая спектральную составляющую с нулевой частотой в начале массива:

module receiver.fft;

private
{
    import std.algorithm;
    import std.complex;
    import std.conv;
    import std.math : log10, PI, cos;
    import std.numeric;
    import std.range;
    import std.string;
}

enum WindowFunction
{
    BLACKMAN,
    BLACKMAN_E,
    BLACKMAN_HARRIS,
    BLACKMAN_NUTALL,
    FLAT_TOP,
    HAMMING,
    HANN,
    NUTALL,
    RECTANGULAR,
    SINE
}

auto blackmanWindow(ref float[] data)
{
    float[] result;

    enum float ALPHA = 0.16f;

    enum float a0 = (1.0f - ALPHA) * 0.5f;
    enum float a1 = 0.5;
    enum float a2 = 0.5 * ALPHA;

    foreach (n, e; data)
    {
        auto tmp = a0 - a1 * cos((2 * PI * n) / (data.length - 1.0)) + a2 * cos((4 * PI * n) / (data.length - 1.0));
        result ~= tmp * e;
    }

    return result;
}

auto blackmanEWindow(ref float[] data)
{
    float[] result;

    enum float a0 = 0.42659;
    enum float a1 = 0.49656;
    enum float a2 = 0.076849;

    foreach (n, e; data)
    {
        auto tmp = a0 - a1 * cos((2 * PI * n) / (data.length - 1.0)) + a2 * cos((4 * PI * n) / (data.length - 1.0));
        result ~= tmp * e;
    }

    return result;
}

auto blackmanHarrisWindow(ref float[] data)
{
    float[] result;

    enum float a0 = 0.35875;
    enum float a1 = 0.48829;
    enum float a2 = 0.14128;
    enum float a3 = 0.01168;

    foreach (n, e; data)
    {
        auto tmp = a0 - a1 * cos((2 * PI * n) / (data.length - 1.0)) + a2 * cos((4 * PI * n) / (data.length - 1.0)) - a3 * cos((6 * PI * n) / (data.length - 1.0));
        result ~= tmp * e;
    }

    return result;
}

auto blackmanNutallWindow(ref float[] data)
{
    float[] result;

    enum float a0 = 0.3635819;
    enum float a1 = 0.4891775;
    enum float a2 = 0.1365995;
    enum float a3 = 0.0106411;

    foreach (n, e; data)
    {
        auto tmp = a0 - a1 * cos((2 * PI * n) / (data.length - 1.0)) + a2 * cos((4 * PI * n) / (data.length - 1.0)) - a3 * cos((6 * PI * n) / (data.length - 1.0));
        result ~= tmp * e;
    }

    return result;
}

auto flatTopWindow(ref float[] data)
{
    float[] result;

    float a0 = 1.0;
    float a1 = 1.93;
    float a2 = 1.29;
    float a3 = 0.388;
    float a4 = 0.028;

    foreach (n, e; data)
    {
        auto tmp = a0 - a1 * cos((2 * PI * n) / (data.length - 1.0)) + a2 * cos((4 * PI * n) / (data.length - 1.0)) - a3 * cos((6 * PI * n) / (data.length - 1.0)) + a4 * cos((8 * PI * n) / (data.length - 1.0));
        result ~= tmp * e;
    }

    return result;
}

auto hammingWindow(ref float[] data)
{
    float[] result;

    foreach (n, e; data)
    {
        auto tmp = 0.53836 - 0.46164 * cos((2 * PI * n) / (data.length - 1.0));
        result ~= tmp * e;
    }

    return result;
}

auto hannWindow(ref float[] data)
{
    float[] result;

    foreach (n, e; data)
    {
        auto tmp = 0.5 * (1.0 - cos((2 * PI * n) / (data.length - 1.0)));
        result ~= tmp * e;
    }

    return result;
}

auto nutallWindow(ref float[] data)
{
    float[] result;

    enum float a0 = 0.355768;
    enum float a1 = 0.487396;
    enum float a2 = 0.144232;
    enum float a3 = 0.012604;

    foreach (n, e; data)
    {
        auto tmp = a0 - a1 * cos((2 * PI * n) / (data.length - 1.0)) + a2 * cos((4 * PI * n) / (data.length - 1.0)) - a3 * cos((6 * PI * n) / (data.length - 1.0));
        result ~= tmp * e;
    }

    return result;
}

auto sineWindow(ref float[] data)
{
    float[] result;

    foreach (n, e; data)
    {
        auto tmp = cos(((PI * n) / data.length) - (0.5 * PI));
        result ~= tmp * e;
    }

    return result;
}

auto applyWindow(ref float[] data, WindowFunction window)
{
    float[] result;

    final switch (window) with (WindowFunction)
    {
        case BLACKMAN:
            result = blackmanWindow(data);
            break;
        case BLACKMAN_E:
            result = blackmanEWindow(data);
            break;
        case BLACKMAN_HARRIS:
            result = blackmanHarrisWindow(data);
            break;
        case BLACKMAN_NUTALL:
            result = blackmanNutallWindow(data);
            break;
        case FLAT_TOP:
            result = flatTopWindow(data);
            break;
        case HAMMING:
            result = hammingWindow(data);
            break;
        case HANN:
            result = hannWindow(data);
            break;
        case NUTALL:
            result = nutallWindow(data);
            break;
        case SINE:
            result = sineWindow(data);
            break;
        case RECTANGULAR:
            result = data;
            break;
    }

    return result;
}

auto makeWindow(string windowName)
{
    WindowFunction window;

    switch (windowName)
    {
        case "Blackman":
            window = WindowFunction.BLACKMAN;
            break;
        case "Blackman(e)":
            window = WindowFunction.BLACKMAN_E;
            break;
        case "Blackman-Nutall":
            window = WindowFunction.BLACKMAN_NUTALL;
            break;
        case "Blackman-Harris":
            window = WindowFunction.BLACKMAN_HARRIS;
            break;
        case "Flat top":
            window = WindowFunction.BLACKMAN_HARRIS;
            break;
        case "Hamming":
            window = WindowFunction.HAMMING;
            break;
        case "Hann":
            window = WindowFunction.HANN;
            break;
        case "Nutall":
            window = WindowFunction.NUTALL;
            break;
        case "Sine":
            window = WindowFunction.SINE;
            break;
        default:
            window = WindowFunction.RECTANGULAR;
            break;
    }

    return window;
}

auto makeFFT(ref float[] data, WindowFunction window)
{
    return data
              .applyWindow(window)
              .fft
              .map!(a => cast(float) (20.0f * log10((a * a).abs)))
              .array;
}

auto centerFFT(float[] data)
{
    auto N = cast(uint) data.length;
    //auto centered = new float[N];

    for (uint i = 0; i < N / 2; i++)
    {
        //centered[i] = data[(N / 2) + i];
        //centered[(N / 2) + i] = data[i];
        data[i] = data[(N / 2) + i];
        data[(N / 2) + i] = data[i];
    }

    //return centered;
    return data[0..$/2];
}

И вот сейчас настало время собрать весь разрозненный микс программных и аппаратных описаний в одну монолитную систему, для чего возвращаемся в файл mainform.d, в котором у нас хранится описание графического интерфейса и некоторых важных функций взаимодействия. Первое, что мы делаем определяем реакцию на кнопку Start: по нажатию на эту кнопку происходит создание объекта устройства USBDevice (это и есть наш COM-порт), затем происходит отключение кнопки Start с одновременным включением кнопки Stop, помимо этого в строку статуса будет выведено сообщение о том, что считывание активно, и сообщение о том, что считывание не удалось, если порт открыть по какой-то причине нельзя. Реакция на нажатие кнопки Stop прямо противоположно описанной реакции нажатия на кнопку Stop.

А теперь самая насыщенная часть нашего проекта — реакция на событие таймера.

И хоть я знаю про то, что смешивать логику проекта с описанием интерфейса нельзя, но тем не менее в этот раз можно с большой натяжкой отступить от этого правила, хоть это и неправильно…

Итак, каждый раз когда срабатывает таймер мы получаем размер области сигнала (а именно длину этой области) и это значение будет использовано в качестве размера блока полученного сигнала. Затем мы создаем перечисление оконной функции с нужной оконной функцией, и считываем первый блок сигнала (т.е получаем первые два байта с устройства и превращаем их в значение сигнала в вольтах), если значение сигнала неравно специальному маркеру -1.0, то добавлем его в массив данных, иначе ничего не делаем. Если длина данных на этот момент равна размеру сигналов, то начинаем обработку данных в цикле: добавляем данные в виджет сигнала, обновляем виджет, выделяем окно соответственно его размеру (размер окна получается из текущего значения слайдера, которое храниться в переменной windowSize), после чего выполняем преобразование Фурье с цетрировкой результата и вводим результат в виджет сонограммы, обновляя его.

Для полноты картины опишем реакцию на изменение в значении слайдера, подставив значение слайдера в выражение, дающее соответствующую степень двойки, сохраним его в переменную windowSize и разместив значение в надпись рядом со слайдером.

Весь код этих процедур выглядит так:

void runStart()
{
    string currentPort = combo0.text!string;

    try
    {
        device = new USBDevice(currentPort);

        startButton.setEnabled(false);
        stopButton.setEnabled(true);
        status0.showMessage("Reading from USB device %s ...".format(currentPort));
    }
    catch (Throwable)
    {
        status0.showMessage("Failed to open USB device %s".format(currentPort));
    }
}

void runStop()
{
    string currentPort = combo0.text!string;

    if (device !is null)
    {
        stopButton.setEnabled(false);
        startButton.setEnabled(true);

        status0.showMessage("Reading from USB device %s stopped".format(currentPort));
        device.close;
    }
}

void runTimer()
{
    auto length = area0.width;
    windowFunction = makeWindow(combo1.text!string);
    float block = device.read;

    if (block != -1)
    {
        data ~= block;
    }

    if (data.length == length)
    {
        for (int i = 0; i < data.length; i++)
        {
            try
            {
                area0.setData(data);
                area0.update;

                auto window = data[i..i + windowSize];
                auto fd = window
                                .makeFFT(windowFunction)
                                .centerFFT;

                area1.setData(fd);
                area1.update;
            }
            catch (Throwable)
            {
                break;
            }
        }

        data = [];
    }
}


void runChange()
{
    windowSize = cast(int) (2 ^^ slider0.value);
    label0.setText("<h3>%d</h3>".format(windowSize));
    label0.update;
}

Испытаем код, скомпилировав с помощью dub и залив предварительно нужный скетч в Arduino (описан выше), а также подключив медный провод к аналоговому выводу A0:

dub run

Результат:

Получилось неплохо, а?!

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