CC BY
124
Сер. 10. 2009. Вып. 1
ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА
УДК 004.383.3(075)
В. С. Тутыгин, А. В. Дебелова
НОВЫЙ АЛГОРИТМ ЧАСТОТНОГО АНАЛИЗА
Исследование ряда объектов методами ядерного магнитного резонанса и ядерного квадрупольного резонанса (ЯМР и ЯКР) приводит к необходимости распознавания слабых эхо-сигналов, лежащих, как правило, значительно (на 20 дБ и более) ниже уровня шумов. Общепринятый подход для распознавания наличия и идентификации эхо-сигнала в зарегистрированном сигнале основан на многократном повторении эксперимента и синхронном накоплении эхо-сигналов с последующим частотным анализом сигнала с помощью быстрого преобразования Фурье (БПФ). Предлагаемый алгоритм частотного анализа зашумленных эхо-сигналов обеспечивает значительно более высокую точность определения частоты.
Задача вычисления точного значения количества периодов и частоты дискретизированного сигнала весьма актуальна при спектральном анализе сигналов спинового и квадрупольного эха импульсных спектрометров ЯМР и ЯКР. Частотный спектр сигнала в этом случае характеризует состав изучаемого вещества, так как каждое вещество (химический элемент) характеризуется строго индивидуальной резонансной частотой. Частотный спектр сигнала нужно определять максимально точно, поскольку частотные линии разных веществ могут находиться весьма близко.
Результат БПФ дискретизированного сигнала определенной частоты - количество периодов сигнала. Если частота отсчетов или интервал дискретности по времени при измерении сигнала известен, то по количеству периодов можно установить и частоту измеряемого сигнала. Точность определения частоты в спектре входного сигнала вполне определенна и зависит от количества периодов p сигнала. Если количество периодов целое, то частота с помощью БПФ находится абсолютно точно (при отсутствии зашумленности сигнала). Если же количество периодов не является целым, то появляется погрешность определения частоты. Максимальное значение погрешности равняется 1/p. В некоторых, практически важных случаях, например, при обработке эхо-сигналов (рис. 1) импульсных спектрометров ЯМР, количество периодов анализируемого сигнала принципиально ограничено величиной около 10 [1]. В этом случае погрешность определения частоты с помощью БПФ достигает 1/10, т. е. 10%.
Влияние шума в регистрируемом сигнале во временной области можно значительно уменьшать за счет многократного повторения эксперимента и синхронного накопления эхо-сигналов (рис. 1).
Однако, как показано в [1], увеличение количества накоплений позволяет улучшать отношение сигнал/шум без искажения формы и уменьшения амплитуды накопленного
Тутыгин Владимир Семенович — доцент кафедры информационных и управляющих систем факультета технической кибернетики Санкт-Петербургского государственного политехнического университета. Количество опубликованных работ: 30. Научное направление: цифровые системы обработки данных и изображений. E-mail: [email protected].
Дебелова Анна Валерьевна — инженер-программист РИВЦ Пулково. Количество опубликованных работ: 3. Научное направление: цифровые системы обработки данных. E-mail: angel [email protected].
© В. С. Тутыгин, А. В. Дебелова, 2009
эхо-сигнала лишь до некоторого предела. После этого предела накопление уже не приносит ощутимого улучшения качества. В частности, эксперименты, проведенные нами на образцах бората железа, показали, что количество накоплений целесообразно брать не более 1000. В данном случае причиной такого ограничения служит то, что резонансная частота бората железа значительно флуктуирует при изменении температуры. Кроме того, при большом времени накопления на эхо-сигнал начинают влиять постепенные изменения параметров приборов, входящих в состав импульсного спектрометра ЯМР. При ограничении времени проведения анализа веществ количество возможных накоплений сигнала должно быть ограничено или вообще должно отсутствовать. Проблема обработки эхо-сигналов в условиях значительных шумов возникает еще и в том случае, когда имеется необходимость исследования веществ в чрезвычайно малых концентрациях. Поэтому задача выполнения частотного анализа зашумленных эхо-сигналов актуальна. В предлагаемом алгоритме упомянутые выше недостатки, присущие БПФ, устранены за счет сочетания БПФ с корреляционным анализом, сплайн-интерполяцией и передискретизацией, что позволило значительно более точно определять значение частоты эхо-сигнала как при целом, так и при нецелом количестве периодов сигнала, при малом количестве периодов сигнала во временном окне и при амплитудной модуляции сигнала.
Рис. 1. Исходный (А) и накопленный (Б) эхо-сигнал в борате железа РеБОз
Конкретная решаемая нами задача частотного анализа рассматривалась применительно к сигналам, получаемым с импульсных спектрометров ЯМР [1] и ЯКР.
Метод ЯКР основан на физическом явлении, типичном для так называемых квад-рупольных ядер в кристаллических структурах. В упорядоченной кристаллической структуре все квадрупольные ядра имеют определенную частоту резонанса, т. е. частоту, на которой происходит резонансное поглощение электромагнитной энергии. Для каждого химического соединения, содержащего такие ядра, существуют одна или несколько характерных резонансных частот, которые определяются структурой данного соединения. По этой причине ЯКР - исключительно избирательный метод для обнаружения тех или иных веществ, содержащих квадрупольные ядра. Особенность ЯКР следующая: в настоящее время исследовано более 10 000 химических соединений, и среди них нет ни одной пары веществ, ЯКР частоты которых совпадали бы. Частоту ЯКР можно считать как бы паспортом химического соединения (вещества) [2, 3].
Предлагаемые новый алгоритм и программа, реализованные в среде MATLAB, решают задачу определения точного значения количества периодов и частоты дискретизированного сигнала как при целом, так и при нецелом количестве периодов и в условиях значительной зашумленности.
Непосредственное использование БПФ не позволяет получить точное значение количества частоты измеряемого сигнала в том случае, когда количество периодов мало, не является целым и/или если обрабатываемый сигнал зашумлен. Алгоритм БПФ, реализованный в среде MATLAB, позволяет установить частоту анализируемого сигнала и в условиях значительных шумов, однако его эффективность значительно снижается, если количество периодов измеряемого сигнала мало и сигнал занимает лишь незначительную часть временной области. Необходимость в изучении именно таких сигналов возникает, в частности, при решении задач обнаружения и идентификации веществ по спектрам ЯКР (рис. 2).
Рис. 2. Временная реализация (А) и спектр (Б) свободной индукции кокаина
В практически важных случаях решение задачи обнаружения конкретного набора веществ и их идентификации по результатам частотного анализа эхо-сигналов спектрометра форма и расположение максимума эхо-сигнала на временной оси по отношению к возбуждающим импульсам и максимуму экспоненциальной огибающей известны. Это дает возможность создать эталонные сигналы, соответствующие ожидаемым для каждого из веществ, и производить идентификацию на основе корреляционного сравнения их с эхо-сигналом. Коэффициент корреляции эхо-сигнала с эталонным будет равен единице, если частоты данных сигналов равны. Однако, так будет только при условии, если эхо-сигнал не зашумлен. Если же он зашумлен, то таким способом определить частоту эхо-сигнала невозможно, так как коэффициент корреляции эхо-сигнала с эталонным будет всегда меньше единицы, причем причиной этого может быть как несогласованность их частот, так и наличие шума в равной степени.
Основа предлагаемого алгоритма заключается в том, что в небольшой окрестности от предполагаемой частоты сигнала (приближенное значение частоты сигнала можно найти с помощью БПФ) вычисляются коэффициенты корреляции с несколькими эталонными сигналами, затем с помощью сплайн-интерполяции и передискретизации
с увеличением количества точек массива строится функция, выражающая зависимость коэффициента корреляции от частоты эталонов. Функция, построенная таким образом, имеет вид параболы с явно выраженным максимумом (рис. 3) как в случае незашум-ленного, так и зашумленного сигнала. При наличии шума форма функции сохраняется, уменьшается лишь абсолютное значение максимума. По положению максимума находится частота эхо-сигнала более точно, чем позволяет сделать БПФ.
Коэф. корреляции ^
1.005
0.995- / Л,
0.99 • J \
0.985- [ \
0.98 - \
0.975- 6 \
0.971_____,_____,______,_____,_____1—
995 1000 1005 1010 1015 1020
Рис. 3. Зависимость коэффициента корреляции от частоты при отсутствии шума (А) и при отношении сигнал/шум 1/3 (Б)
Точное значение частоты равно 1010. Аппроксимирующие кривые построены с помощью функции сплайн-аппроксимации spaps в MATLAB.
Точность определения частоты сигнала тем выше, чем ближе начальное приближение к истинной частоте. Поэтому алгоритм Дебеловой-Тутыгина (ДТ-алгоритм) применяет итерационное вычисление, на каждом этапе итерации в качестве начального приближения используется уточненное значение частоты, полученное на предыдущем этапе. В качестве первого приближения берется частота, определенная с помощью БПФ. Приведем краткое описание алгоритма. Программа, реализующая его, описана в [4].
1. Выполнение БПФ рабочего числового массива исходных данных X[г], i = 0,..., 2n, и получение числового массива частотного спектра Y [j].
2. Нахождение значения jMaKC, соответствующего максимальному значению Y. Найденная величина j соответствует приближенному значению периодов в числовом массиве X [i].
3. Создание 2к + 1 эталонных числовых массивов исходных данных Xs [i], s =
0,1, 2,..., 2k, с числом периодов jMaKC — 1 + s/k.
4. Вычисление коэффициентов корреляции рабочего числового массива со всеми эталонными и формирование числового массива коэффициентов корреляции KK[m], m = 0,1, 2,..., 2k.
5. Выполнение сплайн-интерполяции для массива KK[m] (нахождение непрерывной функциональной зависимости F(m), соответствующей массиву KK[m]).
6. Выполнение передискретизации на основе найденной функциональной зависимости F(m) для массива KK[m] с увеличением количества элементов массива в r раз, т. е. формирование массива KK1 [mi], где m1 = m * r.
7. Нахождение значения m1 макс, соответствующего максимальному значению KK1 (значение m1 макс/г будет представлять вещественное число, отвечающее уточненному
(в общем случае нецелому) значению количества периодов в рабочем числовом массиве исходных данных X [i]).
8. Вычисление разности err = miMaKC/r — jMaKC.
9. Если err < еггдоп, то jMaKC = miMaKC/r и переход к п. 3, иначе переход к п. 10.
10. Вывод найденного точного значения количества периодов miMaKC/r (в общем случае не целого) и частоты f = (m1 MaKC/r)/(dt * 2n) , где dt - шаг дискретности по времени при измерении сигнала X(t).
Точность определения основной частоты при использовании предложенного алгоритма зависит от значений k и r и тем выше, чем они больше, однако, если анализируемый сигнал зашумлен, существенное ее увеличение происходит при росте k и r лишь до некоторого предела. В частности, при соотношении сигнал/шум 1/2... 1/3, k = 3 и r = 10, как показали приведенные исследования, оказывается наилучшим выбором по критерию точность/время анализа.
Количество итераций для расчета частоты с заданной точностью с помощью описанного алгоритма зависит от того, насколько близко к искомой частоте будет находиться начальное приближение. Эффективность же БПФ в качестве средства для определения начального приближения частоты с увеличением уровня зашумленности сигнала снижается, погрешность определения частоты возрастает. Это связано с тем, что числовой массив, представляющий сигнал в частотной области, при зашумленности сигнала во временной области также будет зашумлен и, когда отношение сигнал/шум в частотной области станет больше единицы, правильно установить (хотя бы приблизительно) частоту с помощью БПФ будет невозможно. В связи с этим для определения начального приближения частоты при обработке значительно зашумленных сигналов были предложены FSA-алгоритм и программа [5].
На рис. 4 в качестве примера приведены экранная форма предлагаемой программы и результаты обработки зашумленного сигнала. Амплитуда шума превышает амплитуду сигнала в 2 раза.
Назначение элементов управления и ввода/вывода:
Mode - позволяет выбрать режим работы программы - TestMode или MainMode. В режиме TestMode формирование анализируемого дискретизированного (т. е. представленного в виде числового массива) тестового сигнала производится в программе. В режиме MainMode анализируемый дискретизированный сигнал считывается из файла.
Noise/Signal - окно для ввода отношения шум/сигнал в режиме TestMode.
SignalFrequency - окно для ввода частоты в режиме TestMode. Размерность частоты может быть любой (Гц, кГц, МГц), но она должна совпадать с размерностью частоты отсчетов, задаваемой в окне SamplingFrequency.
SamplingFrequency - окно для ввода частоты отсчетов в режиме TestMode и MainMode. Размерность частоты может быть любой (Гц, кГц, МГц).
FFT - окно для вывода значения частоты, полученной с помощью БПФ. Используется стандартный алгоритм, реализованный в MATLAB.
DT-algorithm - окно для вывода значения частоты, полученной по программе, реализующей ДТ-алгоритм.
При вычислении частоты по массиву данных во временной области, полученному с помощью АЦП с памятью, например LA_h10M6PCI, и занесенному в файл данных, необходимо задать в окне «SamplingFrequency» фактическую частоту отсчетов.
На рис. 5 приведена экспериментально полученная в режиме компьютерного моделирования зависимость погрешности определения частоты зашумленных эхо-сигналов
Рис. 4. Экранная форма предлагаемой программы частотного анализа в среде MATLAB и результаты обработки зашумленного сигнала
в диапазоне частот от 1 до 5 МГц при частоте отсчетов 100 МГц. Генерация нормально распределенного шума выполнена с помощью функции таиви с СКО, равным 0.5, при амплитуде эхо-сигнала, равной 1. Время затухания равнялось 1,3 мкс. Усреднение произведено по 200 реализациям, значение частоты изменялось с шагом 10 кГц.
Погрешность определения частоты после усреднения при применении ДТ-ал-горитма, как следует из рис. 5, в этих условиях не превышает 0,15%. При использовании БПФ погрешность составила 3%. Коэффициент вариации частоты при использовании ДТ-алгоритма составил 1%, при БПФ - 4,5%.
ДТ-алгоритм исследован нами в режиме компьютерного моделирования в среде МЛТЬЛБ с нормальным, равномерным и белым шумом и при обработке эхо-сигналов импульсного спектрометра ЯМР. Для регистрации сигналов применялся быстродействующий АЦП с памятью ЬЛ_н10Ы6РС1 (50 МГц, 2 канала, 256 кбайт памяти). Программа регистрации и обработки данных была реализована в среде LabWiиdows/ СУI 8.0 [6].
Эффективность предлагаемых алгоритма и программы частотного анализа исследовалась в отношении сигналов, подобных получаемым с импульсных спектрометров ЯМР с одной частотной составляющей и экспоненциальной огибающей. В случае, когда сигнал в частотной области имеет несколько хорошо разрешенных максимумов, предлагаемый алгоритм позволяет определить центральную частоту каждого из них.
Погрешность, %
0.2
0.15
0.1
0.05
0
-0.05
-0.1
-0.15
-0.2
Рис. 5. Среднее значение погрешности определения частоты при использовании ДТ-алгоритма, полученное путем компьютерного моделирования.
Частота отсчетов - 100 МГц
Результаты проведенных исследований возможностей нового метода позволили сделать следующие заключения:
1) ДТ-алгоритм дает возможность определять частоту точнее, чем БПФ, даже при большом количестве периодов сигнала во временной области;
2) ДТ-алгоритм позволяет получить точное значение частоты при обработке неза-шумленных сигналов;
3) ДТ-алгоритм по сравнению с БПФ обеспечивает более высокую точность при обработке зашумленных сигналов;
4) ДТ-алгоритм дает возможность определить значительно точнее, чем БПФ, частоту при сверхмалом (от двух) количестве периодов как незашумленного, так и зашумленного сигналов в временной области.
Предлагаемые нами алгоритм и программа частотного анализа спектра, ориентированные на задачи частотного анализа в условиях зашумленности анализируемого сигнала, ограниченности количества периодов сигнала во временной области и амплитудной модуляции, созданы и изучены для решения конкретной задачи, но могут быть применены для частотного анализа других подобных сигналов. Высокая эффективность разработанного нового метода частотного анализа обеспечивается при условии, если количество периодов сигнала во временной области будет не менее двух, а количество отсчетов за период - не менее 20.
Tutygin V. S., Debelova A. V. New algorithm of frequency analysis.
The problem to identify magnetically ordered substances in nuclear magnetic resonance and nuclear quadrupole resonance technologies is usually solved by their echo signal analysis in a
Summary
frequency domain using fast Fourie Transform procedure. In conditions of a low signal to noise ratio this approach brings to significant errors in resonance frequency detection. A new approach is based on echo-signal frequency evaluation using frequency dependence for maximum correlation among an echo-signal specter in frequency domain and several reference specter located in the vicinity of a frequency value, obtained using a usual fast Fourie Transform technique, followed by sharpening this value using a successive approximation principle.
Key words: nuclear magnetic resonance, identification, frequency analyses, fast Fourie Transform, spline-interpolation, rediscretization, cross-correlation.
Литература
1. Тарханов В. И., Тутыгин В. С. Приборный комплекс для поиска и исследования сигналов ЯМР в магнитоупорядоченных веществах // Научное приборостроение. 2003. Т. 13, № 1. С. 58-63.
2. Шелков В. А. Бесконтактный способ выявления взрывчатых и наркотических веществ // http://st.ess.ru/publications/6_2000/shelkov/shelkov.htm.
3. Белый Ю. И., Поцепня О. А., Семин Г. К. и др. Аппаратура для борьбы с терроризмом на основе эффекта ЯКР// http://st.ess.ru/publications/2_2002/semin/semin.htm.
4. Дебелова А. В., Тутыгин В. С. Программа определения точных значений количества периодов и частоты дискретизированного сигнала. Свид. Роспатента № 2007612320 о гос. регистрации программы от 1 июня 2007 г.
5. Тутыгин В. С., Шедов С. В. Программа частотного анализа дискретизированного сигнала. Свид. Роспатента № 2007613363 о гос. регистрации программы от 9 октября 2007 г.
6. Дебелова А. В., Тутыгин В. С. Powerful algorithm of the frequency analysis for LabWindows/CVI 8.0 // http://www.youtube.com/AnnaDeb.
Статья рекомендована к печати проф. Е. И. Веремеем.
Статья принята к печати 7 октября 2008 г.
