Восстановление сигнала по частотному спектру в среде Igor Pro 8

автор: aave

Be the first to comment! Компьютеры и программы Четверг, 28 Ноябрь 2019

Описана методика преобразования спектра сигнала в частотной области в форму сигнала во временной области с помощью обратного быстрого преобразования Фурье и математического пакета Igor Pro 8.

1Описание исходных данных и постановка задачи

Допустим, у нас есть данные о частотном спектре сигнала (трек), полученные с анализатора спектра или любым другим способом. Многие устройства поддерживают язык команд SCPI. В частности, для моего спектроанализатора SCPI команда запроса трека такая:

TRAC1? TRACE1

Необходимо по любому из поддерживаемых интерфейсов (Ethernet, RS-232 и т.д.) послать прибору команду, и если команда подразумевает ответ, прибор пришлёт ответ. Обычно синтаксис команд похож или даже одинаков для устройств аналогичного назначения. Данные приходят в текстовом виде и сохраняются в текстовый или CSV файл. Например, мой спектроанализатор сохраняет трек в таком виде:

27.11.2019 12:32:15	-7.484508514E+001,-7.413254547E+001,-7.459085846E+001,... 
27.11.2019 12:32:20	-7.473252869E+001,-7.492073822E+001,-7.472429657E+001,...

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

Задача состоит в том, чтобы восстановить представление сигнала во временной области по данным спектра. То есть ответить на вопрос: какой сигнал порождает спектр, который мы видим на экране спектроанализатора? Для этого мы будем использовать, конечно же, обратное дискретное преобразование Фурье.

2Восстановление сигнала по спектру в математическом пакете Igor Pro 8

Давайте познакомимся с математическим пакетом Igor Pro 8 фирмы WaveMetrics. В нём имеется богатая функциональность по анализу самых разнообразных данных. Кроме того, эта математическая среда позволяет писать довольно сложные программы (скрипты), которые помогут с обработкой, преобразованием и анализом данных. Подобное мы видели в среде MATLAB (m-файлы).

Запустим программу. Перед нами пустое окно с открытым обозревателем объектов (Data Browser). Здесь будут сохраняться все наши исходные и промежуточные данные, которые к тому же можно структурировать по своему усмотрению. По умолчанию все наши данные будут находиться в корневой ветке root (аналог корня в дереве файлов).

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

Если бы данные располагались в столбец, то их можно было бы просто импортировать: Data Load Waves Load Waves… (можно настроить параметры импорта) или Data Load Waves Load Delimited Text… (для текста с разделителями) или Data Load Waves Load General text… (если в файле содержатся столбцы чисел).

Для написания скрипта откроем окно ввода процедур: Windows Procedure Window (или сочетание Ctrl+M). В этом окне введём код, который здесь представлен. Код обильно снабжён комментариями. Если вкратце, мы открываем файл, считываем его построчно, в каждой строке выделяем числа и заполняем ими массивы (waves).

Код макроса Igor Pro для получения текстовых данных, хранящихся построчно

// Добавляем пункты в меню "Maсros" и назначаем им обработчики:
Menu "Macros"
    "Load lines into waves", /Q, LoadLinesIntoWaves()
    "IFFT", /Q, InverseFFT()
    "FFT", /Q, DirectFFT()
End

// Преобразует файл со строками в массивы.
Function LoadLinesIntoWaves()
   
    variable refNum // хранит ссылку на файл
   
    // выбираемоткрываем окно выбора файла с данными с фильтром по типу файлов
    Open /R /D /F="Все файлы:.*;" refNum 
    if (strlen(S_fileName) == 0)
        return 0 // если ничего не выбрано, выходим из функции
    endif
   
    Open /R refNum as S_fileName // открываем файл и запоминаем ссылку на него
   
    string lineContent // хранит содержимое строки файла
    variable i, j // счётчики цикла
     
    // построчно считываем файл
    string outputWaveName, strValue
    for (i=0; ; i+=1) 
        FReadLine refNum, lineContent // читаем строку из файла и сохраняем в переменную lineContent
        if (strlen(lineContent) == 0)
            return 0 // если больше нет данных, выходим из цикла чтения
        endif

        string firstline = StringFromList(0, lineContent, ",") // первое значение содержит метку времени и измерение
        string timestamp = StringFromList(0, firstline, "\t")
       
        variable nValues = ItemsInList(lineContent, ",") // выделяем из строки данные, разделённые запятыми

        outputWaveName = "myWave" + num2str(i) // сохраняем имя выходного массива
        Make /O /N= (nValues) /D $outputWaveName // создаём выходной массив с заданным именем
        wave output = $outputWaveName  

        // Заполняем массив измерений. Первое значение обрабатываем отдельно,
        // т.к. первое значение содержит метку времени и измерение
        string firstvalue = StringFromList(1, firstline, "\t") 
        output[0] = str2num(firstvalue) // преобразуем строку в число и сохраняем в массив
        
        // Все последующие элементы содержат только данные измерений, проходимся по ним в цикле:
        for (j=1; j<nValues; j+=1)
            strValue = StringFromList(j, lineContent, ",")
            output[j] = str2num(strValue)
        endfor
    endfor
   
    Close refNum
End

Нажмём кнопку Compile, размещённую внизу этого окна. После этого в главном меню программы в разделе Macros появятся три новых пункта: "Load lines into waves", "IFFT" и "FFT". Нажмём первый из них.

Откроется стандартное диалоговое окно выбора файла. Выбираем файл с данными, жмём ОК. Если файл большой, обработка займёт некоторое время. По завершении в просмотрщике объектов появится узел myWave1 (или с другим именем, как вы напишете в скрипте). Если выделить его нажатием, то в нижней части браузера объектов появится предпросмотр данных, которые содержатся в этом узле, в графическом виде. Если сделать на нём двойной щелчок мышью, откроется таблица, где можно посмотреть и отредактировать данные. Если сделать двойной щелчок на графике, откроется окно с полноразмерным графиком и шкалами.

Добавили файл с исходными данными спектра в проект Igor Pro 8

Обратное преобразование Фурье требует, чтобы в данных спектра было число отсчётов, равное степени двойки. В моём случае, как видно из рисунка, 501 отсчёт. Необходимо либо дополнить до 512, либо сократить до 256. Давайте добавим: в данном случае это не сложно, т.к. в конце спектра белый шум. Для этого нажмём на 501-ой строке правой кнопкой мыши и выберем пункт Insert Points…

В Igor Pro нумерация строк идёт с 0, поэтому строке с номером 500 соответствует 501-ая по счёту строка.

Редактирование таблицы данных в Igor Pro 8

Откроется окно Insert Points. Можно добавлять точки по одной, можно сразу указать нужное число точек (Number of points: 11).

Жмём кнопку Do It. В низу таблицы появятся 11 новых строк с нулевыми значениями. Желательно сделать их равными последнему «настоящему» значению, иначе спектр будет искажён. Вообще, в каждом конкретном случае нужно смотреть, какими значениями дополнять. Но тут, как я уже говорил, находятся шумы, и мы просто добавим ещё несколько фиктивных измерений шума с тем же уровнем мощности.

Отлично, можно приступать к восстановлению сигнала во временной области по сигналу в частотной области. В Igor Pro 8 это делается ровно одной строкой кода. Добавим в окне процедур (Procedure Window) ниже функции LoadLinesIntoWaves() следующее.

Макрос Igor Pro для обратого преобразования Фурье

// Обратное БПФ (восстановление сигнала по спектру).
function InverseFFT()
	IFFT /DEST=root:myIfft root:myWave1 
	// root: - это папка по умолчанию; её можно не указывать, если других папок нет;
	// если имеются подпапки, то они должны быть тут добавлены в виде root:subfolder1:subfolder2 и т.п.;
	// если в названии узла имеются пробелы, то оно должно быть заключено в одинарные кавычки, например: root:subfolder1:'my ifft 3'
	// myWave1 - название исходного массива данных.
	// myIfft - имя узла, в котором будет сохранён результат преобразования.
end

Можно ввести код IFFT /DEST=root:myIfft root:myWave1 прямо в консоль (Command Window, которое открывается сочетанием клавиш Ctrl+J) и нажать Enter. Результат будет точно такой же. Но для повторного использования удобно создать функцию и добавить её в меню макросов.

Не забудем скомпилировать код (кнопка Compile в нижней части окна процедур). Выбираем в меню Macros пункт IFFT. В браузере объектов появится новый узел myIfft (или с другим именем, как вы укажете в коде). Если на него установить курсор мыши, снизу отобразятся два графика. Красным цветом показана действительная часть комплексного сигнала, синим – мнимая. Или косинусные и синусные компоненты. Красный график – это и есть сигнал во временной области, который мы хотели восстановить. Именно он даёт исходный частотный спектр.

Результат выполнения обратного преобразования Фурье

Задача решена, но интересно было бы проверить, насколько точно обратное преобразование Фурье восстановило сигнал. Для проверки совершим над полученным сигналом прямое преобразование Фурье. Добавим в Procedure Window следующую функцию:

Макрос Igor Pro для прямого преобразования Фурье

// Прямое БПФ (проверка правильности перевода из спектра в сигнал).
function DirectFFT()
	FFT /DEST=root:myFft root:myIfft 
end

Скомпилируем код, потом выберем в меню Macros пункт FFT. В браузере объектов появится новый узел myFft. Это частотный спектр, полученный из восстановленного сигнала. Он должен совпадать или хотя бы быть похожим на исходный спектр.

Результат выполнения прямого преобразования Фурье над восстановленным сигналом

Как видно, красный график (действительная часть спектра) очень похож на наш исходный спектр. Это значит, что восстановление сигнала по спектру можно считать успешным.

Можно улучшить результат, если сделать следующее. Судя по спектру, он имеет множество мелких возмущений – шумов. Они вносят некторый вклад в суммарную мощность спектра, но не такой большой, как сильные изменения. Можно сгладить спектр перед выполнением обратного преобразования Фурье. Тогда влияние шумов уменьшится. Для этого в браузере объектов выделим узел с нашим исходным сигналом, а затем пойдём в главное меню программы Igor Pro и выберем пункт Analysis Smooth…. В открывшемся диалоговом окне в поле Smoothing введём коэффициент, например, 2. И нажмём кнопку Do it.

Можно поэкспериментировать и подобрать оптимальный коэффициент, при которм шумы будут минимальны, а восстановленный сигнал не будет искажён. В результате этой операции в браузер объектов добавится новый узел со сглаженным спектром: