CC BY
1113MATLAB — новые возможности
в технологии спектроскопии и спектрометрии
Владимир ДЬЯКОНОВ, д. т. н., профессор
Преобразования Фурье сыграли огромную роль в появлении и развитии ряда важных областей науки и техники. Достаточно отметить технику разделения и фильтрации сигналов в радиотехнике и создание современных цифровых анализаторов спектра, в том числе реального времени. Эти очень дорогие приборы аппаратно реализуют быстрое оконное преобразование Фурье и построение трехмерных спектров — спектрограмм [1—3]. Большие и в то же время вполне доступные возможности проведения спектрального анализа сигналов, основанного на этих преобразованиях, открывает матричная система компьютерной математики МА^АВ [4—7], в частности ее новейшие версии МА^АВ 1?2009Ь и МА^АВ 1?2010а. Их возможности в технологии спектроскопии и спектрометрии описаны в этой статье.
Прямое и обратное БПФ в системе МА^АВ
Система компьютерной математики MATLAB, основанная на мощных и современных матричных методах вычислений, является одним из лучших языков программирования высокого уровня для реализации технических вычислений [4]. Она имеет простой и ненавязчивый интерфейс (рис. 1), удобный редактор листингов программ и вычислительное ядро, поддерживающее параллельные вычисления и многоядерные процессоры. Описанные ниже примеры могут быть реализованы как в интерактивном
командном режиме (с набором команд после знака >> в текущей строке командного окна), так и в виде простых программных модулей — m-файлов, редактируемых в текстовом редакторе программ. Фрагменты программ можно вводить в командную строку и в виде копий с текстового процессора Word.
MATLAB содержит функции для выполнения одномерного и двумерного быстрого преобразования Фурье. Для одномерного массива ж с длиной N прямое и обратное преобразования Фурье реализуются по следующим формулам:
Х(к) =J]x(j)e2nmj-l)(k-l\
7=1
Рис. 1. Интерфейс системы компьютерной математики MATLAB
х(к) = —f^X(k)e-2nNIU-m-l\
N k=i
Прямое преобразование Фурье переводит описание сигнала из временной области в частотную, а обратное преобразование Фурье — из частотной области во временную. На этом основаны многочисленные методы фильтрации сигналов.
В описанных далее функциях системы MATLAB реализован особый метод быстрого преобразования Фурье — Fast Fourier Transform (FFT, или БПФ), позволяющий резко уменьшить число арифметических операций в ходе приведенных выше преобразований. Он особенно эффективен, если число обрабатываемых элементов (отсчетов) составляет 2”, где m — целое положительное число.
Для одномерного БПФ в MATLAB используется следующая функция:
• fft(X) — возвращает для вектора Хрезульта-ты его дискретного преобразования Фурье, по возможности используя алгоритм быстрого преобразования Фурье. Если X — матрица, функция fft возвращает преобразование Фурье для каждого столбца матрицы.
• fft(X,n) — возвращает результаты n-точечного преобразования Фурье. Если длина вектора X меньше n, то недостающие элементы заполняются нулями. Если длина X больше n, то лишние элементы удаляются. Когда X — матрица, длина столбцов корректируется аналогично.
• fft(X,[],dim) и fft(X,n,dim) — применяют преобразование Фурье к одной из размерностей массива в зависимости от значения параметра dim.
В качестве иллюстрации типичного применения преобразования Фурье возьмем трехчастотный сигнал на фоне сильного шума, создаваемого генератором случайных чисел (функция randn):
% Задание трехчастотного сигнала с шумом 1[=0:0.0005:1;
х=$т(2*рг*200:Ч)+0.4*5т(2*р1*150:Ч)+0.4*5т(2*р1*250:Ч); y=x+randn(size(t)); рЫ(у(1:2000),Ъ’)
Этот сигнал (рис. 2) содержит 2000 точек и описывает явно синусоидальный сигнал с частотой 200 Гц и два боковых сигнала с частотами 150 и 250 Гц, что соответствует амплитудно-модулированному (АМ) сигналу с частотой модуляции 50 Гц и глубиной модуляции 0,8 (амплитуда боковых частот составляет 0,4 от амплитуды центрального сигнала). После знака % задается неисполняемый комментарий.
На рис. 2 нетрудно заметить, что из него никоим образом не видно, что полезный сигнал — амплитудно-модулированное колебание
и что он содержит синусоидальные компоненты, настолько он забит шумами. И выглядит просто как широкая шумовая дорожка.
Теперь построим график спектральной плотности полученного сигнала с помощью прямого преобразования Фурье. Этот график (рис. 3) в области частот до 300 Гц строится следующим модулем:
% Построение графика спектральной плотности Y=fft(y,1024); Pyy=Y.*conj(Y)/1024; f=2000*(0:150)/1024; plot(f,Pyy(1:151)),grid
Даже беглого взгляда на рис. 3 достаточно, чтобы убедиться в том, что спектр сигнала имеет явный пик на средней (несущей) частоте АМ-сигнала и два боковых пика. Все эти три частотные составляющие сигнала четко выделяются и разделяются на общем шумовом фоне. Из-за наличия шумов пики боковых частот немного разнятся по высоте. Таким образом, данный пример прекрасно иллюстрирует технику обнаружения слабых АМ-сигналов
на фоне шумов, лежащую в основе работы многих радиоприемных устройств.
Для двумерного прямого преобразования Фурье используется функция fft2, а для многомерного — Ап. Функция Y = ДЬМАХ) перегруппировывает выходные массивы функций А и fft2, размещая нулевую частоту в центре спектра, что иногда более удобно. Если X — вектор, то Y — вектор с циклической перестановкой правой и левой половин исходного вектора. Если X — матрица, то Y — матрица, у которой квадранты I и III меняются местами с квадрантами II и IV.
Теперь построим график спектральной плотности мощности (рис. 4) при одномерном преобразовании Фурье для двухчастотного зашумленного сигнала:
% Построение графика спектральной плотности мощности rand (‘state’,0); t=0:0.001:0.512; x=sin (2*pi*50*t)+sin (2*pi*120*t); y=x+2*randn (size (t))+0.3;
Y=fft (y); Pyy=Y.*conj (Y)/512;
f=1000* (0:255)/512; plot (f, Pyy (1:256)), grid
Здесь мы ограничились 29 = 512 отсчетами с тем, чтобы использовать эффективный метод быстрого преобразования Фурье, при котором число отсчетов желательно иметь равным 2”, где N — целое число. Воспользуемся функцией fftshift.
Y=fftshift (Y); Pyy=Y.*conj (Y)/512; plot (Pyy), grid
Полученный при этом график представлен на рис. З. Надо отметить, что этот график дает значения спектральной плотности составляющих спектра не явно от частоты, а как распределение ее значений для элементов вектора Pyy.
Возможно и одномерное обратное преобразование Фурье, реализуемое функцией ifft(F): оно возвращает результат дискретного обратного преобразования Фурье для вектора F. Если F — матрица, то ifft возвращает обратное преобразование Фурье для каждого столбца этой матрицы; ifft(F,n) — возвращает результат n-точечного дискретного обратного преобразования Фурье вектора F:
• ifft(F,[],dim) и y = ifft(X,n,dim) — возвращают результат обратного дискретного преобразования Фурье массива F по строкам или по столбцам в зависимости от значения скаляра dim.
Для любого X результат последовательного выполнения прямого и обратного преобразований Фурье ifft fft (ж)) равен Х с точностью до погрешности округления. Если Х — массив действительных чисел, то ifft (fft (ж)) может иметь мнимые части. Пример (для вектора V):
>> V=[1 1 1 1 О О О О];
>> fft(V); ifft(fft(V)) ans =
1 1 1 1 О О О О
Аналогичные функции есть для двумерного и многомерного случаев.
Математическое моделирование простых сигналов
MATLAB открывает широчайшие возможности в математическом моделировании сигналов. При этом становится предельно ясной математическая сущность сигналов. В MATLAB есть ряд функций, позволяющих задавать сигналы наиболее распространенных форм. Аппаратно эту возможность реализуют генераторы сигналов произвольной формы, но это дорогие приборы [8].
Рассмотрим генерацию импульсных сигналов в MATLAB. Функция y = gmonopuls(t,fc) для заданного вектора отсчетов времени t создает вектор отсчетов y Гауссового двухполярного моноимпульса (рис. 6):
fc = 2E9; fs=100E9; tc = gmonopuls(‘cutoff’, fc); t = -2*tc : 1/fs : 2*tc; y = gmonopuls(t,fc); plot(t,y)
Другая функция у = diric(x,n) служит для задания вектора значений сигнала, представленного функцией Дирихле:
diric(x,n) =
-(л-1)
-I2” ,х=0,±2к,4к...
sin(nx/2)
п sin( X /w)
...else
Параметр п — целое положительное число. Число элементов вектора у равно числу элементов вектора х. Функция diric периодическая, при этом период кратен 2р при нечетных п и 4р при четных п. На рис. 7 показано построение графика сигнала, представленного функцией Дирихле с помощью следующих команд:
х=0:.1:10; subplot(1,1,1); рЫ(х^к(х,4))
Функция gauspuls служит для создания волнового пакета, представляющего собой синусоиду, модулированную по амплитуде
функцией Гаусса. Данная функция может использоваться в нескольких видах. В первом из них уі = gauspuls(t,fc,Ъw[,Ъwr]) она создает вектор отсчетов для моментов времени, заданных в векторе t. Параметр fc задает центральную частоту синусоиды, а Ъw — ширину полосы частот сигнала. По умолчанию заданы ^ = 1000 и Ъw = 0,5. Необязательный параметр Ъwr задает сигнал единичной амплитуды с частотой С и шириной полосы частот bw, причем граница полосы частот задается ослаблением амплитуды на заданное число децибел
Ъwr (по умолчанию----------6 дБ). Этот параметр
должен иметь отрицательное значение.
Следующие две формы функции расширяют ее возможности:
[уі,уд] = gauspuls(...) и [уі^,уе] = gauspuls(...).
В первом случае помимо вектора отсчетов сигнала уі возвращается вектор отсчета дополнительного сигнала, фаза которого сдвинута на 90°. Во втором случае дополнительно к этому возвращается вектор отсчетов огибающей сигнала уе.
Наконец, еще одна форма задания функции:
^ = gauspuls(<cutoff’,fc, bw,bwr,tpe)
Она служит для вычисления времени отсечения С которое определяется по спаду амплитуды до уровня ре дБ (по умолчанию----60 дБ).
Приведем пример применения функции gauspuls:
tc = gauspuls(<cutoff’,З0E3,.б,[],-40); t = -tc : 1E-6 : tc; yi = gauspuls(t,50E3,.75); plot(t,yi)
Нетрудно заметить, что график этой функции (рис. 8) представляет собой волновой пакет и напоминает (чисто внешне) вейвлет. Особенность сигнала этого вида заключается в его временной локализации. Приведенный
Рис. 8. График синусоиды, модулированной функцией Гаусса
Рис. 9. График сигнала, построенного функцией pulstran
пример достаточно прост, но читателю рекомендуется разобраться во всех деталях его синтаксиса.
Функция y = pulstran(t,d, June [,p1,p2,...]) служит для создания отсчетов импульсных сигналов разной формы. Форма задается параметром June, который может иметь значения:
• gauspuls — синусоида, модулированная по закону Гаусса;
• rectpuls — прямоугольный импульс;
• tripuls — треугольный импульс.
Вектор у вычисляется для отсчетов времени, заданных вектором t, по формуле:
у = func(t-d(1)) + func(t-d(2)) + ...
Число импульсов в заданном интервале времени задается длиной вектора d, то есть length (d). Необязательные параметры p1, p2, ... при необходимости позволяют задавать дополнительные параметры обращения к func, например типаfunc(t-d(1),p1,p2,...). При записи функции в виде у = pulstran(t,d,p,[fs]) можно задать частоту дискретизации fs (по умолчанию — 1 Гц).
Ограничимся приведенным ниже примером применения функции pulstran:
t = 0 : .00001 : .005; d = [0 : .001 : .01 ; 0.5.л(0:10)]'; y = pulstran(t,d,@gauspuls,5000,.5); plot(t,y)
График, который построен по этому примеру, показан на рис. 9. Он представляет несколько первых из 10 пакетов (задаются параметром d) гауссовых импульсов, имеющих частоту несущей 5000 Гц и относительную полосу частот 0,5.
Следует обратить особое внимание на возможность этой функции создавать пакеты сигналов и периодические сигналы.
Функция ж = sawtooth (t,[width]) создает вектор пилообразных или треугольных периодических колебаний, уровень которых меняется от -1 до +1 на периоде 2р. Если задано значение параметра width, то им-
пульс на интервале от 0 до 2гех width нарастает в указанных значениях, а на интервале от 2рх width до 2р уменьшается от +1 до -1. При width = 0,5 формируется симметричное треугольное колебание. Приведем пример:
t^^^pi); x1=sawtooth(t,1); x2=sawtooth(t,1/2); plot(t,x1 + 1/2,t,x2-1/2)
Здесь задается построение двух векторов пилообразного х1 и треугольного х2 сигналов (с параметром width, равным 1 и 1/2 соответственно) во временном интервале от 0 до 4р, а затем функция sawtooth строит график этих сигналов (рис. 10). Один из сигналов поднят по оси у на величину 1/2, а другой опущен на эту же величину.
Функция х = square(t,[duty]) генерирует вектор сигнала прямоугольной формы с периодом 2р для моментов времени, имеющихся в векторе t. Положительная полуволна импульсов имеет значение +1, отрицатель-
ная-----1. Параметр duty (по умолчанию ЗО)
задает продолжительность положительной части полуволны колебаний в процентах от периода. К примеру, следующие команды:
t^:.!^; plot(t, square(t, ЗО))
строят график симметричных прямоугольных импульсов, именуемый меандром (рис. 11). Здесь полезно обратить внимание на то, что график построен в пространстве, в котором расположена плоскость графика.
Функция у = tripuls(T[,w[,s]]) служит для создания вектора значений треугольных апериодических импульсов. В форме у=Гр^Т) генерируется одиночный треугольный импульс единичной амплитуды, центрированный относительно Т = 0 и имеющий ширину 1. Параметр w позволяет установить ширину импульса, а параметр s (-1< s <+1) задает асимметрию импульса (по умолчанию s = 0).
Рис. 11. График меандра в пространстве
Рис. 12. График треугольного импульса, полученного с помощью функции tripuls
Следующий вполне очевидный пример:
t=-10:.1:10; plot(t,tripuls(t,З,0.З))
показывает генерацию и построение графика скошенного ^ = 0,5) треугольного импульса с шириной w = 5. Его график в диапазоне времени от -10 до 10 приведен на рис. 12.
Функция sinc задает вектор (или матрицу) сигнала, соответствующего выражению:
• /ч Г1-..*=о
V = вшей) = < .
[вт^О/я*...*Ф 0
Размер вектора (или матрицы) тот же, что у вектора (матрицы) Ь Функция sinc(t) представляет обратное преобразование Фурье для прямоугольного импульса с высотой 1 и шириной 2р:
sinc(0 = J e^dm.
Кроме того, эта функция может использоваться как базисная для восстановления любого сигнала g(t) по его отсчетам при условии, что спектр сигнала ограничен условием -ге<ю<ге:
00
g(.t) = '£g(n)smc(t-n).
Л =-<30
Это положение, вытекающее из известной теоремы Котельникова, иллюстрирует приведенный ниже пример:
randn('state',0); t = (1:13)';
x = [О 1 2 3 2 1 О -1 -2 -3 -2 -1 -О]';
ts = linspace(-З,20,б00)';
y = sinc(ts(:,ones(size(t))) - t(:,ones(size(ts)))')
plot(t,x,'o',ts,y)
Здесь сигнал — одиночный импульс треугольного вида — задан векторами времени (13 отсчетов) и значений параметра (также
13 отсчетов). Функция linspace генерирует 6ОО отсчетов эталонного времени в интервале от -З до 2О. Вектор y задает интерполяционную модель для 13 значений сигнала, описанную выше. Наконец, команда plot строит исходные точки и кривую их интерполяции с помощью функции sinc(t) (рис. 13). Нетрудно заметить, что кривая интерполяции проходит точно только через узловые точки. Этот вид интерполяции широко применяется в современных цифровых осциллографах. При этом не всегда обращают внимание на то, что в промежутках между узлами точное восстановление сигнала не гарантируется.
Функция y = vco(x,fc,fs) создает вектор косинусоидального сигнала с частотной модуляцией (ЧМ). Средняя частота сигнала с единичной амплитудой задается параметром fs. Вектор управляющего воздействия x должен содержать действительные значения воздействия в диапазоне его значений от -1 до +1. При этом отклонение меняется от О до 2/С. Размер вектора y определяется размером вектора x.
Построение спектрограмм сигналов
В форме у = усо (x,[Fmin,Fmax],fs) можно задать изменение частоты от Fmin до Fmax при изменении значений вектора х от -1 до +1. Желательно, чтобы изменение частоты не превышало fs/2. Аргумент ж может быть и матрицей. Тогда описанные правила распространяются на столбцы матрицы, и выход у также будет матрицей. Следующий пример показывает применение данной функции для построения спектрограммы частотно-модулированного сигнала на основе функции specgram:
fs = 1ОО; t = О:.ОО1:2;
x = vco(sawtooth(2*pi*t,0.7З),[0.1 0.4]*fs,fs); specgram(x,512,fs,kaiser(256,5),220)
Эта функция использует метод кратковременного оконного преобразования Фурье, аппаратно реализованный только в дорогих анализаторах спектра реального времени. Вид полученной спектрограммы представлен
на рис. 14. Спектрограмма строится в плоскости время-частота, при этом уровень частотных составляющих сигнала задается цветом. Она отчетливо выявляет характер модуляции по пилообразному закону, тогда как построение графика обычной временной зависимости ЧМ-сигнала этого сделать обычно не позволяет.
Функция у = chirp (t,f0,t1,f1,[‘method’,phi]) формирует выборку (дискретные значения) косинусоидального сигнала с частотой от f0 в начальный момент времени t до f1 в конечный момент времени t1. Звук такого сигнала напоминает визг, откуда и его название. По умолчанию t = 0, f0 = 0 и f1 = 100. Необязательный параметр phi (по умолчанию 0) задает начальную фазу сигнала. Другой необязательный параметр method задает закон изменения частоты. Этот параметр может принимать следующие значения:
• linear — линейный закон изменения частоты fj(t) = f0+at, где a = (f1-f0)/t1;
• quadratic — квадратичный закон изменения частоты f(t) = f0+at2, где a = (f1-f0)/t1;
• logarithmic — логарифмический закон изменения частоты f(t) = f0+10at,
где a = [Iog10(f1-f0)]/t1 и f1>f0.
По умолчанию принято значение method = = linear. Значения параметров по умолчанию используются, если соответствующая переменная отсутствует или задано пустое значение.
Приведем поучительный пример задания косинусоидального сигнала, частота которого меняется по полиномиальному закону. Листинг соответствующего программного модуля выглядит следующим образом:
t=[0 0.5 1.0 1.5 2.0]; % задание вектора времени f=[0 200 100 150 300]; % задание вектора частот p=polyfit(t,f,4); % регрессия полиномом 4-го порядка t=0:0.001:2; % задание вектора времени y=chirp(t,p);%генерация сигнала и построение графиков subplot(211); plot(t,polyval(p,t)); set(gca,'ylim',[0 500]); subplot(212); specgram(y,128,1E3,128,120);
В первых трех строчках модуля задано построение полинома 4-го порядка, описывающего функцию времени, которая используется для модуляции частоты косинусоидального сигнала (следующие две строки). Затем окно графика разбивается на два подокна, и в них строятся графики полинома (рис. 15 сверху) и спектра сигнала (рис. 15 снизу).
Здесь мы видим, что спектрограмма оконного преобразования Фурье прекрасно справляется с задачей идентификации закона частотной модуляции синусоидального сигнала. Это невозможно сделать с помощью осциллографа.
МА^АВ-функции задания окон
Классический спектральный анализ предполагает, что подлежащий анализу сигнал является периодическим и его значения в начале и в конце интервала анализа совпадают. Если это не так, то на спектр сигнала налагаются спектры разрывов, что ведет к появлению лишних составляющих спектра и сильному его искажению.
При обычном оконном спектральном анализе с помощью временных окон из сигнала выделяется область с плавным спадом к нулю у границ интервала анализа. Тем самым разрывы сигнала предотвращаются, хотя сам сигнал несколько искажается. Окна могут быть различного типа и характеризуются графическими зависимостями своих коэффициентов и различными специфическими параметрами. Наиболее широко используются Гауссовы окна, дающие малые искажения спектра сигнала в процессе его ограничения в окнах. Ниже мы рассматриваем окна как зависимость коэффициента передачи окна от номера отсчета к-Щк). Эти коэффициенты являются взвешивающими для корректируемого сигнала.
Пакет расширения Signal Processing Toolbox системы MATLAB имеет ряд функций для задания n-точечных окон. Как правило, они применяются (в том числе в виде вариантов) при выполнении ряда операций спектрального анализа и синтеза. Все функции создают вектор-столбец коэффициентов окна соответствующего типа. Размер его задается параметром n. При n = 1 все функции задания окон возвращают значение 1.
Для окон могут быть построены характеристики амплитудного спектра. Он соответствует частотной характеристике нулевого канала дискретного преобразования Фурье. Для этого может использоваться функция fraqz или просмотр характеристик окна с помощью вьювера vwtool. Фазовые характеристики для всех окон имеют линейный характер и потому особого интереса не представляют.
Вектор w коэффициентов n-точечного окна Бартлетта задается функцией w = bartlett(n). Эти коэффициенты вычисляются по формулам:
• при нечетном n:
w(k) =
2(£-1)/(я-1)К 1<£<(я+1)/2 2—2(£-1)/(и—1)К (п+1)/2<к<п
при четном n:
w(k) =
2(А-1)/(и-1)К \<к<п!2 2(и-£-1)/(и-1)К п/2<к<п-1
Окно Бартлетта подобно треугольному окну, но значение окна Бартлетта при к = 1 и к = 0 равно нулю. Команда:
>> w=bartlett(32);plot(w)
строит окно Бартлетта для n = 32. Оно показано на рис. 1б с открытой позицией меню Tools. Она показывает богатство инструментальных средств окна графики системы MATLAB.
Окно Блэкмана задается функцией w = Ыackman(n,[‘sflag’]). Она возвращает вектор из п коэффициентов данного окна w, вычисляемый по формуле:
w(k) = 0,42-0,5cos
2п
+0,08cos
4л
к-1 и-1 к—І -1
для к = 1, 2, ..., n.
Параметр sflag может иметь следующие значения:
• symmetric — задает симметричное окно (используется по умолчанию);
• periodic — вычисляет окно для (n+1) точки, но возвращает только первые n точек. Команда:
>> w=blackman(32);plot(w)
строит окно Блэкмана для n = 32, показанное на рис. 17.
Функция w = boxcar (n) возвращает n-точечное прямоугольное окно, вычисляемое как w = ones(n,1). Здесь и далее мы не приводим примеры задания окна, поскольку они вполне очевидны и пользователь может составить их по аналогии с приведенными выше примерами.
Окно Чебышева w = chebwin(n,r) задает n-точечный вектор коэффициентов с пульсациями на уровне r (дБ, по умолчанию 100 дБ) в полосе задержания относительно к амплитуде в полосе пропускания. Следующие команды строят график окна Чебышева:
>> w=chebwin(32,120);plot(w)
Окно Чебышева по виду похоже на окно Блэкмана, но имеет немного более острый пик.
№ tat ш «гаї Гак. Lobs
О d Н 4 Г» і Ч'< Й * 1й А * а О L в □
t id*
■ ■
« ■ -
Ы5 ■ ■
Of В7Ї - і -
!■ 1» 15 » В К 3&
Функция w = hamming(n[,‘sflag’]) возвращает вектор w коэффициентов n-точечного окна Хэмминга. Эти коэффициенты вычисляются по формуле:
w(k+\) = 0,54-0,46cos
2я-
п-1
,..к-
= 0.1.
, и-1.
Опция sflag имеет тот же смысл, что и функция задания окна Блэкмана. Читатель может легко построить график этого окна по подобию примеров, приведенных выше.
Функция w = hanning(n[,‘sflag’]) возвращает вектор w коэффициентов п-точечного окна Хэннинга. Эти коэффициенты вычисляются по формуле:
w(k) = 0,5 1—cos
2 71-
п—1
= 1,2, ...,и-1.
...к =
Параметр sflag имеет тот же смысл, что и функция задания окна Блэкмана. График этого окна также несложно построить.
Функция w = kaiser(n,b) задает вектор-столбец п-точечного окна Кайзера. Параметр р задает затухание боковых лепестков окна. Для получения из окна Кайзера фильтра типа КИХ (с конечной импульсной характеристикой) параметр р следует выбирать по формуле:
Р =
0,1102(а- 8,7)... а>50 0,5842(а-21)0>4+ +0,07886(а-21).. ,50>а>21. 0...0<а<21
Следующий пример строит график окна Кайзера (рис. 18).
Рис. 18. Окно Кайзера
>> w=kaiser(32,1);plot(w)
Функция w= Ъгшщ(п) возвращает вектор-столбец коэффициентов п-точечного треугольного окна. При четных п это окно совпадает с окном Бартлетта, за исключением того, что при к = 0 и к = 1 его значение равно 0. При нечетных п коэффициенты треугольного окна вычисляются по формулам:
_ Г2*/(и-1)К 1<£<(п+1)/2 1
[2(т-А:+1)/(и-1)К (я+1)/2<£<и|
Новые функции задания окон
Следующие несколько новых функций служат для задания окон, минимизирующих проявление эффекта Гиббса, обусловленного резким ограничением спектра сигналов при использовании окон:
• w = barthannwin(n) — окно Бартлетта-Хэнна (Вагйей-Напп);
• w = blackmanharris(n) — окно Блэкмана-Харриса (В1асктап-Нагт);
• w = bohmanwin(n) — окно Бохмана (Bohman);
• w = gausswin (n) — окно Гаусса (гауссиана, Gaussian);
• w = gausswin (n,a) — окно Гаусса (Gaussian) с дополнительным параметром a;
• w = nuttallwin(n) — окно Нуталла-Блэк-мана-Харриса (Nuttall’s-Blackman-Harris);
• w = tukeywin (n, a) — окно Тьюкея (Tukey, tapered cosine).
В следующем примере выводятся графики двух окон (blackmanharris и nuttallwin) и строится график их разницы (рис. 19):
N = 64; w = blackmanharris (N); y = nuttallwin (N); subplot(1,2,1); plot(1:N,w,1:N,y,'r--');axis([1 N 0 1]); title('Comparison of 64-pt windows'); legend('Blackman-harris', 'Nuttall'); subplot(1,2,2); plot(y-w);
title('Difference between Blackman-harris and Nuttall windows')
Новая обобщенная функция задания окон w = window (fhandle,n) возвращает n-точечное
окно любого типа, которое задается параме-тром^а^іе (дескриптором), содержащим символ @ и имя окна:
@barthannwin @hamming @bartlett @hann
@blackman @kaiser @blackmanharris @nuttallwin
@bohmanwin @rectwin @chebwin @triang
@gausswin @tukeywin
В приведенном примере строятся графики для трех окон, построенные функцией window (рис. 20).
Итак, MATLAB дает обширные возможности в задании и исследовании окон.
Построение графиков амплитудного спектра окон
В анализаторах спектра окно определяет вид пиков спектра, наблюдаемого на экране анализатора. Для построения графика амплитудного спектра можно использовать функцию freq, что иллюстрирует следующий пример (рис. 21) для окна Хэмминга:
Рис. 21. График амплитудного спектра для окна Хэмминга
w = hamming(20); w=w/sum(w); [h,f]=freqz(w,1,[],20); plot(f, 20*log10(abs(h))); grid on
Полезно обратить внимание на то, что относительные амплитуды нормированы и спектр построен в логарифмических единицах (децибелах). Заменив функцию задания окна в первой строке приведенного фрагмента кода, можно построить графики амплитудного спектра и для других видов окон.
Применение вьювера окон VWTool
В реализации пакета Signal Processing Toolbox 6.0/6.1 появилась новая функция — wvtool(w1[,w2,w3,...wn]), позволяющая в окне вьювера просматривать графики окон w1, w2, ... wn и их амплитудных спектров. Так, команда:
>> wvtool(hann(32),hamming(32),hanning(32))
выводит графики трех окон (Хэнна, Хэмминга и Хэннинга) и их амплитудных спектров, представленные на рис. 22.
В справке по пакету все функции окон имеют графическое представление, полученное с помощью вьювера окон.
Спектральная оценка зашумленных сигналов
В радиотехнике особый интерес представляет спектральная оценка сильно зашумленных сигналов. Для таких сигналов используются два подхода:
• непараметрический — использующий только информацию, извлеченную из сигнала (реализован в методах периодограмм и Уэлча);
• параметрический — предполагающий наличие некоторой статистической модели сигналов, параметры которой подлежат определению (реализован в других методах этого раздела).
Спектр дискретного сигнала является периодическим, и прямое дискретное преобразование Фурье — ДПФ (Discrete Fourier Transform, DFT) определяется выражением:
ІУ — 1
Х(п) = У\(£)ехр
. 2ппк
~J
N
N
і N12-1 _
x(t)=— Y, д«)єхр
Для предотвращения растекания (размазывания) спектра дискретных сигналов часто используются описанные выше окна. Для их применения достаточно в формуле прямого ДПФ под знаком суммы ввести еще один множитель — W(k). Соответственно, обратное дискретное преобразование Фурье задается выражением:
1 N_1 ■ ( Jvnt-
x(k)=—TJx(n)exp\j
Для получения полосы частот сигнала от 0 до р/Т приходится смещать нумерацию отсчетов. При нечетном числе отсчетов суммирование ведется при п, меняющемся от -(”-1)/2 до (”-1)/2. Коэффициенты Х(п) с отрицательными номерами вычисляются из соотношения симметрии.
Частотным спектром случайного процесса является преобразование Фурье от корреляционной функции случайного процесса Кх:
ВД = Ё Д#)ехр(-./а>И’).
Методы спектрального оценивания зашумленных сигналов
функции имеют внутреннее обращение к функции рШ(^Рхх) и строят графики зависимости спектральной плотности мощности (СПМ) от частоты, используя векторы частот f и СПМ на них Рхх.
Непараметрические методы — периодограмм и Уэлча
Периодограммой называют оценку СПМ одной реализации случайного процесса:
2
Г(ю) = — %
N-1
^ х{к) ехр(-у'со&Г)
*=0
При использовании весовых функций (окон) периодограмма вычисляется как:
ДПФ имеет следующие свойства:
• является линейным преобразованием;
• дает задержку на один такт при умножении спектра на множитель exp (-j 2ren/N);
• обладает симметрией, то есть:
X(N-n) = X(-n) =x\n);
• применимо к произведению последовательностей отсчетов одинаковой длины;
• допускает операцию круговой свертки (перемножения спектров) двух последовательностей.
ДПФ легко обеспечивает восстановление непрерывных периодических сигналов с ограниченным спектром. Для этого достаточно номер отсчета k заменить на нормированное время t/T. Тогда формула восстановления при четном числе отсчетов будет иметь вид:
MATLAB с пакетом расширения Signal Processing Toolbox имеет с десяток методов спектрального оценивания зашумленных Щсо) = сигналов. Основные данные по ним представлены в таблице. Все описанные ниже
Таблица. Основные методы спектрального оценивания сигналов
^х(к)Щк) ехр(-jakT)
к=0
fu
Y}w(k)\
NT
Метод Функция Достоинства Недостатки
Периодограмм, непараметрический periodogram Хорошее разрешение по частоте Изрезанность при длинных сигналах. Флюктуации при большом числе отсчетов
Уэлча, непараметрический pwelch Хорошее сглаживание сигналов Возможная потеря острых пиков спектра
Берга, на основе авторегрессии pburg Высокое разрешение и стабильность при анализе коротких сигналов. Хорошее устранение шума. Использование окон (весовых функций) и частично перекрывающихся сегментов сигнала. Дисперсия оценки спектра мощности меньше, чем при методе Бартлетта Начальные фазы синусоид сильно влияют на положение спектральных пиков; появляется смещение спектральных пиков при анализе суммы синусоид с шумом; при большом порядке модели может наблюдаться расщепление спектральных пиков
Юла-Уокера, на основе авторегрессии Корректный учет краевых отсчетов сигнала и минимизация погрешности линейного предсказания. Хорошие результаты при анализе спектров длинных последовательностей при стабильном рассчитанном формирующем фильтре Плохие результаты при анализе коротких сигналов. При анализе зашумленных синусоидальных сигналов возникает смещение центральных пиков
Ковариационный pcov Высокая разрешающая способность при спектральной оценке коротких сигналов. Минимизация погрешностей прямого предсказания. Хорошая оценка сигналов в виде суммы чистых синусоид Рассчитанный формирующий фильтр может оказаться нестабильным. При анализе суммы синусоид с шумом возможно смещение спектральных пиков
Модифицированный ковариационный pmcov Высокое разрешение при анализе коротких сигналов, отсутствие расщепления пиков Зависимость положения спектральных линий от начальных фаз синусоид
Многооконный, на основе последовательностей Слепиана pmtm Высокое разрешение Работа только с вещественными данными
Для построения периодограмм как зависимостей спектральной мощности в дБ/Вт от частоты служит функция periodogram:
Частота может задаваться как угловая w или в герцах f Параметр range может иметь два значения:
• twosided — вычисляет двухполосную спектральную мощность в частотном диапазоне [0,fs). При задании вместо этой опции пустого вектора [] частотный диапазон определяется как [0,1). Если не используется спецификация f, то этот диапазон выбирается равным [0,2).
• onesided — вычисляется однополосная спектральная мощность в диапазоне частот, за-
данном для вещественных компонент вектора х. Для х с вещественными компонентами эта опция используется по умолчанию.
Для построения периодограмм (рис. 23) используется функция periodogram, что иллюстрирует следующий пример:
Fs = 1000; t = 0:№:.3;
х = sm(2*pi*t*200)+cos(2*pi*t*350)+0.01*randn(size(t));
periodogram(x,[],’onesided’,512,Fs)
Главный недостаток периодограммы — ее сильная изрезанность.
Функция у = goertzel(х,) для вектора х возвращает результат его дискретного Фурье-преобразования, используя одноименный с ней алгоритм второго порядка (goertzel). Если х — матрица, то преобразование выполняется для каждого столбца. Вектор должен содержать целые числа от 1 до Ы, где N — значение первого размера для матрицы х, который больше 1.
Метод Уэлча (рис. 24) является улучшенным вариантом метода периодограмм. Он реализуется следующей функцией:
[Рхх^] = pwelch(x) [Рхх^] = (х,пшп) ...
Здесь целочисленное положительное значение параметра nwin задает длину окна Хэмминга. Если nwin задан как двухэлементный вектор, он задает размеры прямоугольного окна. Для расширения знаний об этой функции рассмотрим пример:
Fs = 1000; t = 0:1^:.3; х = cos(2*pi*t*200) + randn(size(t)); pwelch(x,33,32,[],Fs,’twosided’)
Дополнительное описание алгоритма реализации этого метода можно найти в справочной системе МА^АВ.
Параметрические методы спектрального оценивания
Из параметрических методов спектрального оценивания в первую очередь стоит упомянуть метод Берга. Для иллюстрации оценки СПМ методом Берга воспользуемся следующим примером. Вначале построим АЧХ и ФЧХ (рис. 25) автокорреляционной AR-системы:
а = [1 -2.2137 2.9403 -2.1697 0.9606]; %Коэффициенты AR-фильтра freqz(1,a) % АЧХ и ФЧХ AR-фильтра
Теперь построим АЧХ для сигнала, полученного из белого шума фильтраций с помощью AR-системы 4-го порядка методом Берга (рис. 26):
randn('state',1); % Задание белого шума х = filter(1,a,randn(256,1)); % Выход AR-системы pbuгg(x,4) % Оценка СКМ 4-го порядка методом Берга
С работой ковариационного метода можно ознакомиться, выполнив приведенный ниже пример:
U J Ц * I h Ч £ 9 »! - а О в D
* f и і'1 3 5 -Id я
і І і і
\ і
у \ Г
С 3 : _
ЗІ 3 7 9 і Пі Qi 0F §? І.І ^ПСМЧГу ;.| ГЯІЧжфаІ! В V
№ І ■ I I 1 1 !
]
1 "1ІЄ ■1Н |HRn.«n.nIni.mn.n
Г“— : Г ^
I г
3 И 91 б 2- 4 4 ЭЛ М IF Cl U
Рис. 25. АЧХ и ФЧХ AR-фильтра
[Pxx,w] = periodogram(x)
[Pxx,w] = periodogram(x,window) [Pxx,w] = periodogram(x,window,nfft) [Pxx,f] = periodogram(x,window,nfft,fs) [Pxx,...] = periodogram(x,..., ‘range’) periodogram(...)
a = [1 -2.2137 2.9403 -2.1б97 0.9б0б]; freqz(1,a) title('AR System Frequency Response') randn('state',1); x = filter(1,a,randn(256,1)); pmcov(x,4)
Для реализации многооконного метода служит функция pmtm. Параметр nw этой функции задает спектрально-временное разрешение. По умолчанию nw = 4, рекомендуется выбирать его значения равными 2, 5/2, 3, 7/2. Число последовательности Слепиана равно
2nw-1. Функция имеет параметр method, позволяющий задать метод вычисления СПМ:
• adapt — адаптивный нелинейный алгоритм Томсона комбинации индивидуальных оценок (используется по умолчанию);
• unity — линейная комбинация индивидуальных оценок с весами, равными 1;
• eigen — линейная комбинация индивидуальных оценок с весами, задаваемыми собственными значениями.
Приведем пример оценки СПМ синусоиды с частотой 3000 Гц, на которую наложен небольшой шум (рис. 27):
U J . і
Глек fed л
Ч \ О З
/-г - а о Q и a
G,4 ЛЧ-Гу ЕЯНЧИ
!
/ І \ Ї / / 1 \ І / т”"1Ymmm— V
5 14 \/ ] S І і r-^: і E \ і V ]
1 : * : 1 \ І 1 ] Е \ * \і і
* 1 І і І
! “ 1
51 91 її 41 * ]!■ &.Е І F 01 S1 1
Ьят-ден г№іжу ;-і гай^жприї
Рис. 26. АЧХ сигнала с белым шумом, построенная с помощью функции pburg
randn('state',0); fs = 1000; t = 0:1/fs:0.3; x = cos(2*pi*t*300) + 0.1*randn(size(t)); [Pxx,Pxxc,f] = pmtm(x,3.5,512,fs,0.99); hpsd = dspdata.psd([Pxx Pxxc],'Fs',fs); plot(hpsd)
Функция pmusic реализует алгоритм MUSIC (Multiple Signal Classification). Функция pmusic делит автокорреляционную матрицу на два подпространства (две матрицы) для сигнала и шума, в основе чего лежит анализ собственных значений матриц. Функция записывается в виде:
[S,w] = pmusic(x,p)[S,w] = pmusic(...,nfft)
[S,f] = pmusic(x,p,nfft,fs)) ...
Здесь вместо параметра Рхх используется параметр 5, что связано с тем, что данный алгоритм вычисляет псевдоспектр сигнала. Параметр х может быть как вектором, так и матрицей. Если х — вектор, то в нем задаются отсчеты одного сигнала, а если матри-
ца — то отсчеты ряда сигналов по столбцам. Параметр corr означает, что в векторе ж задается корреляционная матрица. Скалярный параметр p определяет размерность подпространства сигнала. Он может быть также двухэлементным вектором [p thresh], причем параметр thresh задает порог выделения подпространств сигнала и шума.
В наиболее полной форме:
[Sf] = pmusic(x,p,nfft,fs,nwin,noverlap)
функции pmusic возможно задание числа отсчетов nfft для БПФ и частоты квантования fs. По умолчанию nfft = 25б и fs = 1 Гц. При вещественном x оценка СПМ ведется в диапазоне частот от 0 до fs/2, а если x содержит комплексные значения, то в диапазоне частот от 0 до fs. Параметр windows может быть скаляром или вектором. Соответственно, он
задает либо длину окна, либо коэффициенты для произвольно задаваемого окна. Наконец, скалярный параметр noverlap задает количество отсчетов для перекрытия окон.
С другими нюансами применения функции pmusic можно познакомиться, используя команду help music, а также справочную систему пакета. В качестве примера оценим СКМ двухчастотного сигнала с шумом:
randn('state',1); n=0:99;
s=exp(i*pi/3*n)+ exp(i*pi/5*n)+randn(1,100);
X=corrmtx(s,12,'mod'); pmusic(X,3,'whole')
Построенный по этому примеру спектр показан на рис. 28.
Для вычисления векторов круговых частот и мощностей (в дБ) служит функция:
[w,pow] = гоо^ш^хф) [f,pow] = гооШшк(...^)...
*г- Ь«« Ти* СпАщд «Ии- "*
dju-і k -a
! Iі I | 1 і ! і і і і -Н——f--" і ! і і J і 1“^ ■■■■■■ : 1 [
4 : 3 j і
г II 1 Г1 ! 1
і Г| і І ІІ.І4 і! ! 3 "
і /,! S £ I І ї tr ■ і Vi : ■ :. ; 1 і .-W rt. ..11-. J
%: їм і»
ЙН ПЙ X* 'r*S***C[f{4f|
:-л іс< i'.i r.M
Рис. 29. Пример оценки СПМ методами Уэлча и Берга для зашумленного двухчастотного сигнала
Ограничимся приведением вполне очевидного примера использования данной функции:
randn('state',1); n=0:99;
s = exp(i*pi*n)+2*exp(i*pi/2*n)+exp(i*pi/3*n)+randn(1,100); X=corrmtx(s,12,'mod'); [W,P] = rootmusic(X,3)
>>rmusic W =
1.5721
1.0492
3.1414
P =
4.3117
1.4105
1.0084
Сравнение спектральных оценок разными методами
Представляет интерес сравнение различных методов спектральной оценки при использовании реальных сигналов. В качестве примера (файл cs) сравним метод Берга с методом Уэлча на сигнале, представленном двумя синусоидами с амплитудами 1 и 2 и частотами 160 и 140 Гц (рис. 29):
randn('state',0); fs = 1000; t = (0:fs)/fs;
A = [1 2];f = [160;140];
xn = A*sin(2*pi*f*t) + 0.1*randn(size(t));
[P1,f] = pwelch(xn,hamming(256),128,1024,fs);
[P2,f] = pburg(xn,14,1024,fs); plot(f,10*log10(P1),':',f,10*log10(P2)); grid ylabel('PSD Estimates (dB/Hz)'); xlabel('Frequency (Hz)'); legend('Welch','Burg')
По аналогии с этим примером читатель может самостоятельно сравнить между собой и другие методы спектрального оценивания. Множество примеров этого можно найти в справке по пакету и в разделе демонстрационных примеров по нему. Далее в статье можно найти интересные примеры спектрального оценивания (анализа) реальных осциллограмм, получаемых от цифровых электронных осциллографов.
MATLAB-функция построения спектрограмм — specgram
Для визуализации короткооконного БПФ служит функция построения спектрограмм
specgram. Это очень мощное и современное средство визуализации спектра, аппаратно реализованное в анализаторах спектра реального времени фирмы Tektronix [2, 3]. Спектрограмма представляется зависимостью амплитуды спектральных составляющих БПФ, вычисляемого в перемещающемся окне, от момента времени, задающего положение окна. Фактически спектрограмма строится в плоскости частота-время, а амплитуда каждой спектральной составляющей определяет цвет построения каждой точки спектрограммы. При построении спектрограммы используется функциональная окраска, иногда с применением шкалы цветов.
Данная функция имеет ряд форм записи. В простейшем случае B = specgram (а) вычисляется спектрограмма сигнала, с отчетами в векторе ж. При этом ряд параметров используется по умолчанию:
• nfft = min (256, length (a)); fs = 2;
• window — периодическое окно Хэннинга с длиной nfft и numoverlap = length (window)/2. В других формах записи могут задаваться
различные входные параметры и определяться дополнительные выходные параметры:
B = specgram(a,nfft)
[B,f] = specgram(a,nfft,fs)
[B,f,t] = specgram(a,nfft,fs)
B = specgram(a,nfft,fs,window[,numoverlap])
B = specgram(a,f,fs,window,numoverlap)
Если какой-то из параметров не задается, то использование пустого списка [] задает его значение по умолчанию. Поскольку назначение всех входных параметров уже не раз обсуждалось, отметим, что наряду с амплитудами спектральных составляющих B может возвращаться вектор частот БПФ f и вектор времен t. Длина вектора t равна числу столбцов матрицы B. Параметр numoverlap задает число отсчетов, на которое происходит перекрытие блоков.
При ж с комплексными компонентами матрица B будет также содержать комплексные
компоненты с числом строк nttf. Каждый столбец соответствует определенному моменту времени, пропорциональному номеру столбца.
Функция specgram (...) строит спектрограмму в текущем окне, используя функцию:
imagesc(t,f,20*log10(abs(b))), axis xy, colormap(jet)
Рассмотрим пару примеров построения спектрограмм. В одном из них строится спектрограмма звуковых колебаний из тестового файла mtlb:
load mtlb; specgram(mtlb,512,Fs,kaiser(500,5),475) title('Spectogramm for audio wave')
Этот простой пример считывает с жесткого диска звуковой файл mtlb и затем строит его спектрограмму, используя для этого окна Кайзера. Полученная спектрограмма представлена на рис. 30.
Алгоритм вычисления спектрограмм содержит три характерных шага:
1. Разбивка x на перекрывающиеся блоки, на каждый из которых накладывается окно.
2. Выполнение nttf-точечного БПФ для соответствующего отрезка времени, что создает соответствующий столбец матрицы B, после чего окно перемещается на число точек, равное (length (window) - numoverlap). Если число точек БПФ превышает количество отсчетов в окне, то перед выполнением БПФ блок дополняется нулями.
3. При вещественных компонентах x спектрограмма строится для положительных частот, и матрица B содержит при nfft четном nfft/2+1 строк, а при nfft нечетном — (nfft+1)/2 строк и k = fix((n-numoverlap)/ /(length(window)-numoverlap)) столбцов. Подобные спектрограммы широко применяются в электроакустике, поскольку позволяют создавать красочный графический образ звуковых колебаний, в котором
О ? " * *& t 1 17 11 111 II t І Jf
Рис. 31. Спектрограмма сигнала vcosig
я* U J4- Лтт\ ҐиЛ Е*іі л- ■*#«■ н* '
ц .о»с/-;а □ е по
ПКГГ44-РІ сити* ирсс <■-*■->: і-гц
з н:« .ом іж> 4'Яї vs «сси іти: -се о жв і&ак
МК^ТИОЧІ АІУІ
Рис. 32. Построение четырех спектров масс-спектрометра
опытный взгляд может подметить множество тонких особенностей анализируемых звуков. Как шутку, в которой много правды, приведем в заключение спектрограмму колебания из файла vcosig (рис. 31), входящего в набор демонстрационных файлов системы МА^АВ. Эта спектрограмма построена с помощью команд:
Обработка в МА^АВ спектрометрических данных в пакете биоинформатики
Одни из самых сложных спектров — спектры, получаемые масс-спектрометрами при исследовании состава веществ и биологических объектов. Их обработка является одним из основных методов контроля состава веществ в спектроскопии и спектрометрии,
а порою и обнаружения новых веществ в микроскопических количествах. Ныне это крайне актуально при борьбе с терроризмом, когда нужно обнаруживать и распознавать малые дозы взрывчатых веществ и вещественные доказательства преступлений. Средства для работы с микромассивами веществ и обработки масс-спектрометрической информации сосредоточены в новом пакете расширения системы MATLAB по биоинформатике (Bioinformatics Toolbox) [9].
Массовые данные спектроскопии и спектрометрии обычно сохраняются в файлах текстового формата CVS, хранящих данные о массе/нагрузке (M/Z) и значениях интенсивности, соответствующих этим отношениям. Для загрузки данных используют одну из функций MATLAB — ввода/вывода: importdata, csvread или textread. Можно использовать также функцию jcampread, чтобы загрузить специальные JCAMP-DX отформатированные файлы, и xlsread, чтобы загру-
зить файлы электронных таблиц в формате рабочих книг Excel.
Следующий пример показывает считывание тестового файла масс-спектрометра и построение графиков четырех спектров (рис. 32), данные которых имеются в четырех файлах формата CVS:
sample = importdata('mspec01.csv');
MZ = sample.data(:,1); Y = sample.data(:,2); files = {'mspec01.csv','mspec02.csv','mspec03.csv','mspec04.csv'}; for i = 1:4 D = csvread(files{i},1); Y(:,i) = D(:,2); end plot(MZ,Y); axis([0 10000 -20 105])
xlabel('Mасса/Заряд (M/Z)');ylabel('Огносительная интенсивность') ШкСЗагрузка и построение спектров масс-спектрометрв')
Спектры легко различаются по цвету. По пикам их кривых можно судить о составе веществ, для которых они получены. Для этого обычно требуется обработка спектров, например, пере-выборка отчетов, коррекция базовой линии, поиск и сортировка пиков спектральных линий и т. д. Все эти возможности реализованы функциями пакета расширения Bioinformatics Toolbox. Следующий пример строит спектр с базовой линией (красной на рис. 33):
YB = msbackadj(MZ,Y,'WIND0WSIZE',500,'QUANTILE',0.20, 'SHOWPLOT',2);
Калибровка масс-спектрометра ведет к изменениям отношения между наблюдаемым M/Z вектором и истинным временем пролета ионов. Погрешности могут накапливаться при повторных экспериментах. Для устранения этого можно использовать функцию msalign, что позволит стандартизировать значения M/Z. В следующем примере задан массив спектрограммы и весов, используемых для осуществления алгоритма сглаживания:
PP = [3991.4 4598 79б4 91б0]; W = [бО 100 бО 100];
Теперь построим карту спектра (по сути, спектрограмму) до выравнивания, используя команды:
load vcosig; specgram(vcosig,[],Fs)
«• ЕЯ. Pt- М "иА Стй-.я ■** ■ r*m ’
^ j и - ч ч ■■. о Ф *> / - a □ £ n □
\—y їж '.вса ic<n rtaa esc vat 'Mg
RteWOwga W7)
Рис. 34. Карта спектра (спектрограмма) после выравнивания
МюаЭДи АВД
Рис. 35. Спектр после обработки
msheatmap(MZ,YB,'MARKERS',P,'LIMIT',[3000 10000]) й^єСПослє выравнивания (Alignment)')
Следующие команды обеспечивают построение карты спектра после процедуры выравнивания (рис. 34):
Для уменьшения погрешностей при обработке спектра используется операция его нормализации. Интенсивность каждого сигнала приводится к некоторому значению, например 100. При этом можно исключить некоторые области спектра. Для нормализации используется функция т$пот.
Основную информацию в спектре веществ несут пики: каждый пик соответствует тому или иному химическому элементу, а высота пика несет информацию о его количестве в веществе. Поэтому конечной целью спектрального анализа является выявление пиков в спектрах веществ и их обработка.
Для выявления пиков спектрограмм можно использовать функцию численного дифференцирования МА^АВ diff. Пик обнаруживается, если численная производная спектра вначале положительна, а затем отрицательна. При
этом порою надо не учитывать пики слишком малой высоты и пики в низкочастотной области спектра (рис. 35). Реализация алгоритма поиска пиков с применением функции find представлена следующим примером:
YA = msalign(MZ,YB,P,'WEIGHTS',W); % выравнивание YN1 = msnorm(MZ,YA,'QUANTILE',1,'LIMITS',...
[1000 inf], 'MAX',100); % нормализация
YN2 = msnorm(MZ,YA,'LIMITS',[1000 inf],'MAX',100); % нормализация
YS = mssgolay(MZ,YN1,'SPAN',35,'SHOWPLOT',3); % снижение шума slopeSign = diff(YS(:,1))>0; % поиск пиков по смене знака slopeSignChange = diff(slopeSign)<0; % производной h = find(slopeSignChange) + 1; % нахождение всех пиков h(MZ(h)<1500)=[]; % устранение низкочастотных пиков h(YS(h,1)<5)=[]; % устранение малых по высоте пиков plot(MZ,YS(:,1),'-',MZ(h),YS(h,1),'ro') % Построение спектра axis([0 10000 -5 100]) % задание масштаба по осям
xlabel('Масса/Заряд (M/Z)');
ylabel('Относительная интенсивность');
ШкСАвтоматическая селекция пиков спектра')
Пакет имеет основанный на графическом интерфейсе пользователя GUI вьювер спектров масс-спектрометров. Вот пример его вывода функцией msviewer:
msviewer (MZ,YS,'MARKERS',MZ(h),'GROUP',1:4)
Спектры во вьювер могут загружаться из рабочего пространства MATLAB или из файлов формата CVS.
В последнее время помимо Фурье-спектро-грамм сигналов стали использоваться вейвлет-спектрограммы, которые имеют ряд преимуществ при анализе тонких особенностей сигналов. Техника вейвлет-преобразования сигналов и ее инструментальные средства (включая вейвлет-спектрограммы) описана в [5, 10]. Там же детально описан пакет расширения системы MATLAB по вейвлетам — Wavelet Toolbox. Применение техники вейвлет-преобразований (и Фурье-преобразований) при анализе реальных осциллограмм от современных цифровых осциллографов описано в [7]. ■
Литература
1. Афонский А. А., Дьяконов В. П. Цифровые анализаторы спектра, сигналов и логики. М.: СОЛОН-Пресс, 2009.
2. Анализаторы спектра реального времени. Tektronix — www.tektonix.com/rsa
3. Дьяконов В. П. Современные цифровые анализаторы спектра // Компоненты и технологии. 2010. № 5.
4. Дьяконов В. П. MATLAB R2006/2007/2008 + Simulink 5/6.7. Основы применения. 2-е изд., дополн. и перераб. М.: СОЛОН-Пресс, 2008.
5. Дьяконов В. П. Современные методы Фурье-и вейвлет-анализа и синтеза сигналов // Контрольно-измерительные приборы и системы. 2009. № 2.
6. Дьяконов В. П. Компьютерная математика в измерительной технике // Контрольноизмерительные приборы и системы. 2009. № 5.
7. Дьяконов В. П. MATLAB — новые возможности в технологии осциллографии // Компоненты и технологии. 2009. № 10.
8. Дьяконов В. П. Генерация и генераторы сигналов. М.: ДМК-Пресс, 2008.
9. Дьяконов В. П., Круглов В. В. MATLAB 6.5 SP1/ 7/7 SP1/7SP2 + Simulink 5/6. Инструменты искусственного интеллекта и биоинформатики. М.: СОЛОН-Пресс, 2006.
10. Дьяконов В. П. Вейвлеты. От теории к практике.
2-е изд., перераб. и дополн. М.: СОЛОН-Р, 2004.
msheatmap(MZ,YB,'MARKERS',P,'LIMIT',[3000 10000]) title('Before Alignment')
