MATLAB — новые возможности в технологии спектроскопии и спектрометрии

Прямое и обратное БПФ в системе МАTLАВ

Система компьютерной математики МАTLАВ, основанная на мощных и современных матричных методах вычислений, является одним из лучших языков программирования высокого уровня для реализации технических вычислений [4]. Она имеет простой и ненавязчивый интерфейс (рис. 1), удобный редактор листингов программ и вычислительное ядро, поддерживающее параллельные вычисления и многоядерные процессоры. Описанные ниже примеры могут быть реализованы как в интерактивном командном режиме (с набором команд после знака >> в текущей строке командного окна), так и в виде простых программных модулей — от-файлов, редактируемых в текстовом редакторе программ. Фрагменты программ можно вводить в командную строку и в виде копий с текстового процессора Word.

Рис. 1. Интерфейс системы компьютерной математики MATLAB

MATLAB содержит функции для выполнения одномерного и двумерного быстрого преобразования Фурье. Для одномерного массива ж с длиной N прямое и обратное преобразования Фурье реализуются по следующим формулам:

Прямое преобразование Фурье переводит описание сигнала из временной области в частотную, а обратное преобразование Фурье — из частотной области во временную. На этом основаны многочисленные методы фильтрации сигналов.

В описанных далее функциях системы MATLAB реализован особый метод быстрого преобразования Фурье — Fast Fourier Transform (FFT, или БПФ), позволяющий резко уменьшить число арифметических операций в ходе приведенных выше преобразований. Он особенно эффективен, если число обрабатываемых элементов (отсчетов) составляет 2^m, где m — целое положительное число.

Для одномерного БПФ в MATLAB используется следующая функция:

• fft(X) — возвращает для вектора Х-результаты его дискретного преобразования Фурье, по возможности используя алгоритм быстрого преобразования Фурье. Если X — матрица, функция fft возвращает преобразование Фурье для каждого столбца матрицы.
• fft(X,n) — возвращает результаты n-точечного преобразования Фурье. Если длина вектора X меньше n, то недостающие элементы заполняются нулями. Если длина X больше n, то лишние элементы удаляются. Когда X — матрица, длина столбцов корректируется аналогично.
• fft(X,[],dim) и fft(X,n,dim) — применяют преобразование Фурье к одной из размерностей массива в зависимости от значения параметра dim.

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

% Задание трехчастотного сигнала с шумом
t=0:0.0005:1;
x=sin(2*pi*200*t)+0.4*sin(2*pi*150*t)+0.4*sin(2*pi*250*t);
y=x+randn(size(t)); plot(y(1:2000),'b')

Этот сигнал (рис. 2) содержит 2000 точек и описывает явно синусоидальный сигнал с частотой 200 Гц и два боковых сигнала с частотами 150 и 250 Гц, что соответствует амплитудно-модулированному (АМ) сигналу с частотой модуляции 50 Гц и глубиной модуляции 0,8 (амплитуда боковых частот составляет 0,4 от амплитуды центрального сигнала). После знака % задается неисполняемый комментарий.

Рис. 2. Форма зашумленного сигнала

На рис. 2 нетрудно заметить, что из него никоим образом не видно, что полезный сигнал — амплитудно-модулированное колебание и что он содержит синусоидальные компоненты, настолько он забит шумами. И выглядит просто как широкая шумовая дорожка.

Теперь построим график спектральной плотности полученного сигнала с помощью прямого преобразования Фурье. Этот график (рис. 3) в области частот до 300 Гц строится следующим модулем:

% Построение графика спектральной плотности
Y=fft(y,1024); Pyy=Y.*conj(Y)/1024;
=2000*(0:150)/1024; plot(f,Pyy(1:151)),grid

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

Рис. 3. График спектральной плотности сигнала, приведенного на рис. 2

Для двумерного прямого преобразования Фурье используется функция fft2, а для многомерного — Ап. Функция Y = ДЬМАХ) перегруппировывает выходные массивы функций А и размещая нулевую частоту в центре спектра, что иногда более удобно. Если 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

Рис. 4. График спектральной плотности сигнала после одномерного преобразования Фурье

Здесь мы ограничились 2⁹ = 512 отсчетами с тем, чтобы использовать эффективный метод быстрого преобразования Фурье, при котором число отсчетов желательно иметь равным 2^м, где N — целое число. Воспользуемся функцией fftshift:

Y=fftshift (Y); Pyy=Y.*conj (Y)/512; plot (Pyy), grid

Полученный при этом график представлен на рис. 5. Надо отметить, что этот график дает значения спектральной плотности составляющих спектра не явно от частоты, а как распределение ее значений для элементов вектора Pyy.

Рис. 5. График спектральной плотности того же сигнала после применения функции fftshift

Возможно и одномерное обратное преобразование Фурье, реализуемое функцией ifft(F): оно возвращает результат дискретного обратного преобразования Фурье для вектора F. Если F — матрица, то ifft возвращает обратное преобразование Фурье для каждого столбца этой матрицы; fft(F,n) — возвращает результат n-точечного дискретного обратного преобразования Фурье вектора F:

• ifft(F,[],dim) и y = ifft(X,n,dim) — возвращают результат обратного дискретного преобразования Фурье массива F по строкам или по столбцам в зависимости от значения скаляра dim.

Для любого X результат последовательного выполнения прямого и обратного преобразований Фурье fftfft(x)) равен Х с точностью до погрешности округления. Если Х — массив действительных чисел, то ifftfft(x)) может иметь мнимые части. Пример (для вектора V):

>> V=[1 1 1 1 0 0 0 0];
>> fft(V); ifft(fft(V))
ans =
1 1 1 1 0 0 0 0

Аналогичные функции есть для двумерного и многомерного случаев.

Математическое моделирование простых сигналов

МА^АВ открывает широчайшие возможности в математическом моделировании сигналов. При этом становится предельно ясной математическая сущность сигналов. В МАTLАВ есть ряд функций, позволяющих задавать сигналы наиболее распространенных форм. Аппаратно эту возможность реализуют генераторы сигналов произвольной формы, но это дорогие приборы [8].

Рассмотрим генерацию импульсных сигналов в МАTLАВ. Функция у = gmonopuls(t,fc) для заданного вектора отсчетов времени t создает вектор отсчетов у Гауссового двухпо-лярного моноимпульса (рис. 6):

fc = 2E9; fs=100E9; tc = gmonopuls('cutoff', fc);
t = -2*tc : 1/fs : 2*tc; y = gmonopuls(t,fc); plot(t,y)

Рис. 6. График Гауссова моноимпульса

Другая функция у = diric(x,n) служит для задания вектора значений сигнала, представленного функцией Дирихле:

Параметр n — целое положительное число. Число элементов вектора у равно числу элементов вектора х. Функция diric периодическая, при этом период кратен 2р при нечетных п и 4р при четных п. На рис. 7 показано построение графика сигнала, представленного функцией Дирихле с помощью следующих команд:

x=0:.1:10; subplot(1,1,1); plot(x,diric(x,4))

Рис. 7. График сигнала, представленного функцией Дирихле при n = 4

Функция gauspuls служит для создания волнового пакета, представляющего собой синусоиду, модулированную по амплитуде функцией Гаусса. Данная функция может использоваться в нескольких видах. В первом из них yi = gauspuls(t,fc,bw[,bwr]) она создает вектор отсчетов для моментов времени, заданных в векторе t. Параметр С задает центральную частоту синусоиды, а bw — ширину полосы частот сигнала. По умолчанию заданы fc = 1000 и bw = 0,5. Необязательный параметр Ъwr задает сигнал единичной амплитуды с частотой fc и шириной полосы частот bw, причем граница полосы частот задается ослаблением амплитуды на заданное число децибел bwr (по умолчанию—6 дБ). Этот параметр должен иметь отрицательное значение.

Следующие две формы функции расширяют ее возможности:

[yi,yq] = gauspuls(...) и [yi,yq,ye] = gauspuls(...).

В первом случае помимо вектора отсчетов сигнала уі возвращается вектор отсчета дополнительного сигнала, фаза которого сдвинута на 90°. Во втором случае дополнительно к этому возвращается вектор отсчетов огибающей сигнала уе.

Наконец, еще одна форма задания функции:

tc = gauspuls('cutoff',fc, bw,bwr,tpe)

Она служит для вычисления времени отсечения С которое определяется по спаду амплитуды до уровня tpe дБ (по умолчанию — –60 дБ).

Приведем пример применения функции gauspuls:

tc = gauspuls('cutoff',50E3,.6,[],-40);
t = -tc : 1E-6 : tc; yi = gauspuls(t,50E3,.75); plot(t,yi)

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

Функция y = pulstran (t,d, ‘func’ [,p1,p2,…]) служит для создания отсчетов импульсных сигналов разной формы. Форма задается параметром func, который может иметь значения:

• gauspuls — синусоида, модулированная по закону Гаусса;
• ectpuls — прямоугольный импульс;
• tripuls — треугольный импульс.

Вектор y вычисляется для отсчетов времени, заданных вектором t, по формуле:

y = func(t-d(1)) + func(t-d(2)) + …

Число импульсов в заданном интервале времени задается длиной вектора d, то есть length (d). Необязательные параметры p1, p2, … при необходимости позволяют задавать дополнительные параметры обращения к func, например типа func(t-d(1),p1,p2,…). При записи функции в виде y = 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.

Рис. 9. График сигнала, построенного функцией pulstran

Следует обратить особое внимание на возможность этой функции создавать пакеты сигналов и периодические сигналы.

Функция x = sawtooth(t,[width]) создает вектор пилообразных или треугольных периодических колебаний, уровень которых меняется от –1 до +1 на периоде 2p. Если задано значение параметра width, то импульс на интервале от 0 до 2p*width нарастает в указанных значениях, а на интервале от 2p*width до 2p уменьшается от +1 до –1. При width = 0,5 формируется симметричное треугольное колебание. Приведем пример:

t=(0:.05:4*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 до 4p, а затем функция sawtooth строит график этих сигналов (рис. 10). Один из сигналов поднят по оси y на величину 1/2, а другой опущен на эту же величину.

Рис. 10. Графики сигналов, построенных функцией sawtooth

Функция х = square(t,[duty]) генерирует вектор сигнала прямоугольной формы с периодом 2р для моментов времени, имеющихся в векторе t. Положительная полуволна импульсов имеет значение +1, отрицательная -1. Параметр duty (по умолчанию 50) задает продолжительность положительной части полуволны колебаний в процентах от периода. К примеру, следующие команды:

t=0:.1:20; plot(t, square(t, 50))

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

Рис. 11. График меандра в пространстве

Функция у = tripuls(T[,w[,s]]) служит для создания вектора значений треугольных апериодических импульсов. В форме у=Мр^Т) генерируется одиночный треугольный импульс единичной амплитуды, центрированный относительно Т = 0 и имеющий ширину 1. Параметр w позволяет установить ширину импульса, а параметр s (-1< s <+1) задает асимметрию импульса (по умолчанию s = 0).

Следующий вполне очевидный пример:

t=-10:.1:10; plot(t,tripuls(t,5,0.5))

показывает генерацию и построение графика скошенного (s = 0,5) треугольного импульса с шириной w = 5. Его график в диапазоне времени от -10 до 10 приведен на рис. 12.

Рис. 12. График треугольного импульса, полученного с помощью функции tripuls

Функция sinc задает вектор (или матрицу) сигнала, соответствующего выражению:

Размер вектора (или матрицы) тот же, что у вектора (матрицы) t. Функция sinc(t) представляет обратное преобразование Фурье для прямоугольного импульса с высотой 1 и шириной 2р:

Кроме того, эта функция может использоваться как базисная для восстановления любого сигнала g(t) по его отсчетам при условии, что спектр сигнала ограничен условием –p<w<p:

Это положение, вытекающее из известной теоремы Котельникова, иллюстрирует приведенный ниже пример:

randn('state',0); t = (1:13)';
x = [0 1 2 3 2 1 0 -1 -2 -3 -2 -1 -0]';
ts = linspace(-5,20,600)';
y = sinc(ts(:,ones(size(t))) - t(:,ones(size(ts)))')*x;
plot(t,x,'o',ts,y)

Здесь сигнал — одиночный импульс треугольного вида — задан векторами времени (13 отсчетов) и значений параметра (также 13 отсчетов). Функция linspace генерирует 600 отсчетов эталонного времени в интервале от -5 до 20. Вектор y задает интерполяционную модель для 13 значений сигнала, описанную выше. Наконец, команда plot строит исходные точки и кривую их интерполяции с помощью функции sinc(t) (рис. 13). Нетрудно заметить, что кривая интерполяции проходит точно только через узловые точки. Этот вид интерполяции широко применяется в современных цифровых осциллографах. При этом не всегда обращают внимание на то, что в промежутках между узлами точное восстановление сигнала не гарантируется.

Рис. 13. Интерполяция сигнала, представленного набором узловых точек

Функция y = vco(x,fc,fs) создает вектор косинусоидального сигнала с частотной модуляцией (ЧМ). Средняя частота сигнала с единичной амплитудой задается параметром fs. Вектор управляющего воздействия x должен содержать действительные значения воздействия в диапазоне его значений от -1 до +1. При этом отклонение меняется от 0 до 2fc. Размер вектора y определяется размером вектора x.

Построение спектрограмм сигналов

В форме y = vco (x,[Fmin,Fmax],fs) можно задать изменение частоты от Fmin до Fmax при изменении значений вектора x от -1 до +1. Желательно, чтобы изменение частоты не превышало fs/2. Аргумент ж может быть и матрицей. Тогда описанные правила распространяются на столбцы матрицы, и выход y также будет матрицей. Следующий пример показывает применение данной функции для построения спектрограммы частотно-модулированного сигнала на основе функции specgram:

fs = 100; t = 0:.001:2;
x = vco(sawtooth(2*pi*t,0.75),[0.1 0.4]*fs,fs);
specgram(x,512,fs,kaiser(256,5),220)

Эта функция использует метод кратковременного оконного преобразования Фурье, аппаратно реализованный только в дорогих анализаторах спектра реального времени. Вид полученной спектрограммы представлен на рис. 14. Спектрограмма строится в плоскости время-частота, при этом уровень частотных составляющих сигнала задается цветом. Она отчетливо выявляет характер модуляции по пилообразному закону, тогда как построение графика обычной временной зависимости ЧМ-сигнала этого сделать обычно не позволяет.

Рис. 14. Спектрограмма частотно-модулированного сигнала

Функция y = chirp(t,f0,t1,f1,[‘method’,phi]) формирует выборку (дискретные значения) косинусоидального сигнала с частотой от f0 в начальный момент времени t до f1 в конечный момент времени t1. Звук такого сигнала напоминает визг, откуда и его название. По умолчанию t = 0, f0 = 0 и f1 = 100. Необязательный параметр phi (по умолчанию 0) задает начальную фазу сигнала. Другой необязательный параметр method задает закон изменения частоты. Этот параметр может принимать следующие значения:

• linear — линейный закон изменения частоты f((t) = f0+at, где a = (f1-f0)/t1;
• quadratic — квадратичный закон изменения частоты f(t) = f0+at2, где a = (f1-f0)/t1;
• logarithmic — логарифмический закон изменения частоты fj(t) = f0+10at, где a = [log10(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 снизу).

Рис. 15. График модулирующей полиномиальной функции и спектрограмма сигнала, модулированного по заданному этой функцией закону

Здесь мы видим, что спектрограмма оконного преобразования Фурье прекрасно справляется с задачей идентификации закона частотной модуляции синусоидального сигнала. Это невозможно сделать с помощью осциллографа.

MATLAB-функции задания окон

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

При обычном оконном спектральном анализе с помощью временных окон из сигнала выделяется область с плавным спадом к нулю у границ интервала анализа. Тем самым разрывы сигнала предотвращаются, хотя сам сигнал несколько искажается. Окна могут быть различного типа и характеризуются графическими зависимостями своих коэффициентов и различными специфическими параметрами. Наиболее широко используются Гауссовы окна, дающие малые искажения спектра сигнала в процессе его ограничения в окнах. Ниже мы рассматриваем окна как зависимость коэффициента передачи окна от номера отсчета k-W(k). Эти коэффициенты являются взвешивающими для корректируемого сигнала.

Пакет расширения Signal Processing Toolbox системы MATLAB имеет ряд функций для задания n-точечных окон. Как правило, они применяются (в том числе в виде вариантов) при выполнении ряда операций спектрального анализа и синтеза. Все функции создают вектор-столбец коэффициентов окна соответствующего типа. Размер его задается параметром п. При n = 1 все функции задания окон возвращают значение 1.

Для окон могут быть построены характеристики амплитудного спектра. Он соответствует частотной характеристике нулевого канала дискретного преобразования Фурье. Для этого может использоваться функция fraqz или просмотр характеристик окна с помощью вьювера vwtool. Фазовые характеристики для всех окон имеют линейный характер и потому особого интереса не представляют.

Вектор w коэффициентов n-точечного окна Бартлетта задается функцией w = bartlett(n). Эти коэффициенты вычисляются по формулам:

• при нечетном n:

• при четном n:

Окно Бартлетта подобно треугольному окну, но значение окна Бартлетта при k = 1 и k = 0 равно нулю. Команда:

>> w=bartlett(32);plot(w)

строит окно Бартлетта для n = 32. Оно показано на рис. 16 с открытой позицией меню Tools. Она показывает богатство инструментальных средств окна графики системы MATLAB.

Рис. 16. Окно Бартлетта

Окно Блэкмана задается функцией w = blackman(n,[‘sflag’]). Она возвращает вектор из п коэффициентов данного окна и, вычисляемый по формуле:

для k = 1, 2, …, n.

Параметр sflag может иметь следующие значения:

• symmetric — задает симметричное окно (используется по умолчанию);
• periodic — вычисляет окно для (n+1) точки, но возвращает только первые n точек.

Команда:

>> w=blackman(32);plot(w)

строит окно Блэкмана для n = 32, показанное на рис. 17.

Рис. 17. Окно Блэкмана для n = 32

Функция w = boxcar (n) возвращает n-точечное прямоугольное окно, вычисляемое как w = ones(n,1). Здесь и далее мы не приводим примеры задания окна, поскольку они вполне очевидны и пользователь может составить их по аналогии с приведенными выше примерами.

Окно Чебышева w = chebwin(n,r) задает n-точечный вектор коэффициентов с пульсациями на уровне r (дБ, по умолчанию 100 дБ) в полосе задержания относительно к амплитуде в полосе пропускания. Следующие команды строят график окна Чебышева:

>> w=chebwin(32,120);plot(w)

Окно Чебышева по виду похоже на окно Блэкмана, но имеет немного более острый пик.

Функция w = hamming(n[,’sflag’]) возвращает вектор w коэффициентов n-точечного окна Хэмминга. Эти коэффициенты вычисляются по формуле:

Опция sflag имеет тот же смысл, что и функция задания окна Блэкмана. Читатель может легко построить график этого окна по подобию примеров, приведенных выше.

Функция и = hanning(n[,’sflag’]) возвращает вектор и коэффициентов п-точечного окна Хэннинга. Эти коэффициенты вычисляются по формуле:

Параметр sflag имеет тот же смысл, что и функция задания окна Блэкмана. График этого окна также несложно построить.

Функция w = kaiser(n,b) задает вектор-столбец п-точечного окна Кайзера. Параметр b задает затухание боковых лепестков окна. Для получения из окна Кайзера фильтра типа КИХ (с конечной импульсной характеристикой) параметр р следует выбирать по формуле:

Следующий пример строит график окна Кайзера (рис. 18).

>> w=kaiser(32,1);plot(w)

Рис. 18. Окно Кайзера

Функция w= triang(n) возвращает вектор-столбец коэффициентов n-точечного треугольного окна. При четных n это окно совпадает с окном Бартлетта, за исключением того, что при k = 0 и k = 1 его значение равно 0. При нечетных n коэффициенты треугольного окна вычисляются по формулам:

новые функции задания окон

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

• w = barthannwin (n) — окно Бартлетта-Хэнна (Bartlett-Hann);
• w = blackmanharris (n) — окно Блэкмана-Харриса (Blackman-Harris);
• w = bohmanwin(n) — окно Бохмана (Bohman);
• w = gausswin (n) — окно Гаусса (гауссиана, Gaussian);
• w = gausswin (n,a) — окно Гаусса (Gaussian) с дополнительным параметром a;
• w = nuttallwin (n) — окно Нуталла-Блэк-мана-Харриса (NuttalTs-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')

Рис. 19. Графики двух окон (слева) и их разность (справа)

Новая обобщенная функция задания окон w = window (fhandle,n) возвращает n-точечное окно любого типа, которое задается параметром fhandle (дескриптором), содержащим символ @ и имя окна:

@barthannwin @hamming @bartlett       @hann
@blackman    @kaiser  @blackmanharris @nuttallwin
@bohmanwin   @rectwin @chebwin        @triang 
@gausswin    @tukeywin

В приведенном примере строятся графики для трех окон, построенные функцией window (рис. 20).

Рис. 20. Графики трех окон, построенные функцией window

Итак, MATLAB дает обширные возможности в задании и исследовании окон.

Построение графиков амплитудного спектра окон

В анализаторах спектра окно определяет вид пиков спектра, наблюдаемого на экране анализатора. Для построения графика амплитудного спектра можно использовать функцию freq, что иллюстрирует следующий пример (рис. 21) для окна Хэмминга:

w = hamming(20); w=w/sum(w); [h,f]=freqz(w,1,[],20);
plot(f, 20*log10(abs(h))); grid on

Рис. 21. График амплитудного спектра для окна Хэмминга

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

Применение вьювера окон VWTool

В реализации пакета Signal Processing Toolbox 6.0/6.1 появилась новая функция — wvtool(w1[,w2,w3,…wn]), позволяющая в окне вьювера просматривать графики окон w1, w2, … wn и их амплитудных спектров. Так, команда:

>> wvtool(hann(32),hamming(32),hanning(32))

выводит графики трех окон (Хэнна, Хэм-минга и Хэннинга) и их амплитудных спектров, представленные на рис. 22.

Рис. 22. Окно вьювера окон

В справке по пакету все функции окон имеют графическое представление, полученное с помощью вьювера окон.

Спектральная оценка зашумленных сигналов

В радиотехнике особый интерес представляет спектральная оценка сильно зашумленных сигналов. Для таких сигналов используются два подхода:

• непараметрический — использующий только информацию, извлеченную из сигнала (реализован в методах периодограмм и Уэлча);
• параметрический — предполагающий наличие некоторой статистической модели сигналов, параметры которой подлежат определению (реализован в других методах этого раздела).

Спектр дискретного сигнала является периодическим, и прямое дискретное преобразование Фурье — ДПФ (Discrete Fourier Transform, DFT) определяется выражением:

Для предотвращения растекания (размазывания) спектра дискретных сигналов часто используются описанные выше окна. Для их применения достаточно в формуле прямого ДПФ под знаком суммы ввести еще один множитель — Щ(К). Соответственно, обратное дискретное преобразование Фурье задается выражением:

• • ДПФ имеет следующие свойства: является линейным преобразованием;
  • дает задержку на один такт при умножении спектра на множитель exp (–j2pn/N);

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

ДПФ легко обеспечивает восстановление непрерывных периодических сигналов с ограниченным спектром. Для этого достаточно номер отсчета k заменить на нормированное время t/T. Тогда формула восстановления при четном числе отсчетов будет иметь вид:

Для получения полосы частот сигнала от 0 до p/Т приходится смещать нумерацию отсчетов. При нечетном числе отсчетов суммирование ведется при п, меняющемся от –(N–1)/2 до (N–1)/2. Коэффициенты X(n) с отрицательными номерами вычисляются из соотношения симметрии.

Частотным спектром случайного процесса является преобразование Фурье от корреляционной функции случайного процесса R_x:

Методы спектрального оценивания зашумленных сигналов

MATLAB с пакетом расширения Signal Processing Toolbox имеет с десяток методов спектрального оценивания зашумленных сигналов. Основные данные по ним представлены в таблице. Все описанные ниже функции имеют внутреннее обращение к функции plot(f,Pxx) и строят графики зависимости спектральной плотности мощности (СПМ) от частоты, используя векторы частот f и СПМ на них Рхх.

Непараметрические методы — периодограмм и Уэлча

Периодограммой называют оценку СПМ одной реализации случайного процесса:

При использовании весовых функций (окон) периодограмма вычисляется как:

Для построения периодограмм как зависимостей спектральной мощности в дБ/Вт от частоты служит функция periodogram:

[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(…)

Частота может задаваться как угловая w или в герцах f Параметр range может иметь два значения:

• twosided — вычисляет двухполосную спектральную мощность в частотном диапазоне [0,fs). При задании вместо этой опции пустого вектора [] частотный диапазон определяется как [0,1). Если не используется спецификация f, то этот диапазон выбирается равным [0,2).
• onesided — вычисляется однополосная спектральная мощность в диапазоне частот, заданном для вещественных компонент вектора ж. Для ж с вещественными компонентами эта опция используется по умолчанию.

Для построения периодограмм (рис. 23) используется функция periodogram, что иллюстрирует следующий пример:

Fs = 1000; t = 0:1/Fs:.3;
x = sin(2*pi*t*200)+cos(2*pi*t*350)+0.01*randn(size(t));
periodogram(x,[],'onesided',512,Fs)

Рис. 23. Построение спектра двухчастотного сигнала методом периодограмм

Главный недостаток периодограммы — ее сильная изрезанность.

Функция у = goertzel(x,i) для вектора ж возвращает результат его дискретного Фурье-преобразования, используя одноименный с ней алгоритм второго порядка (goertzel). Если ж — матрица, то преобразование выполняется для каждого столбца. Вектор должен содержать целые числа от 1 до N где N — значение первого размера для матрицы ж, который больше 1.

Метод Уэлча (рис. 24) является улучшенным вариантом метода периодограмм. Он реализуется следующей функцией:

[Pxx,w] = pwelch(x) [Pxx,w] = (x,nwin) …

Рис. 24. Пример оценки СПМ методом Уэлча

Здесь целочисленное положительное значение параметра nwin задает длину окна Хэмминга. Если nwin задан как двухэлементный вектор, он задает размеры прямоугольного окна. Для расширения знаний об этой функции рассмотрим пример:

Fs = 1000; t = 0:1/Fs:.3; x = cos(2*pi*t*200) + randn(size(t));
pwelch(x,33,32,[],Fs,'twosided')

Дополнительное описание алгоритма реализации этого метода можно найти в справочной системе МАTLАВ.

Параметрические методы спектрального оценивания

Из параметрических методов спектрального оценивания в первую очередь стоит упомянуть метод Берга. Для иллюстрации оценки СПМ методом Берга воспользуемся следующим примером. Вначале построим АЧХ и ФЧХ (рис. 25) автокорреляционной AR-системы:

a = [1 -2.2137 2.9403 -2.1697 0.9606]; %Коэффициенты AR-фильтра
freqz(1,a) % АЧХ и ФЧХ AR-фильтра

Рис. 25. АЧХ и ФЧХ AR-фильтра

Теперь построим АЧХ для сигнала, полученного из белого шума фильтраций с помощью AR-системы 4-го порядка методом Берга (рис. 26):

randn('state',1); % Задание белого шума
x = filter(1,a,randn(256,1)); % Выход AR-системы
pburg(x,4) % Оценка СКМ 4-го порядка методом Берга

Рис. 26. АЧХ сигнала с белым шумом, построенная с помощью функции pburg

С работой ковариационного метода можно ознакомиться, выполнив приведенный ниже пример:

a = [1 -2.2137 2.9403 -2.1697 0.9606]; 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):

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)

Рис. 27. Пример оценки СПМ с помощью функции pmtm

Функция pmusic реализует алгоритм MUSIC (Multiple Signal Classification). Функция pmusic делит автокорреляционную матрицу на два подпространства (две матрицы) для сигнала и шума, в основе чего лежит анализ собственных значений матриц. Функция записывается в виде:

[S,w] = pmusic(x,p)[S,w] = pmusic(...,nfft)
[S,f] = pmusic(x,p,nfft,fs)) …

Здесь вместо параметра Pxx используется параметр 5, что связано с тем, что данный алгоритм вычисляет псевдоспектр сигнала. Параметр х может быть как вектором, так и матрицей. Если х — вектор, то в нем задаются отсчеты одного сигнала, а если матрица — то отсчеты ряда сигналов по столбцам. Параметр corr означает, что в векторе ж задается корреляционная матрица. Скалярный параметр p определяет размерность подпространства сигнала. Он может быть также двухэлементным вектором [p thresh], причем параметр thresh задает порог выделения подпространств сигнала и шума.

В наиболее полной форме:

[S,f] = pmusic(x,p,nfft,fs,nwin,noverlap)

функции pmusic возможно задание числа отсчетов nfft для БПФ и частоты квантования fs. По умолчанию nfft = 256 и 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.

Рис. 28. Спектр двухчастотного сигнала с шумом, построенный функцией pmusic

Для вычисления векторов круговых частот и мощностей (в дБ) служит функция:

[w,pow] = rootmusic(x,p) [f,pow] = rootmusic(...,fs)…

Ограничимся приведением вполне очевидного примера использования данной функции:

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')

Рис. 29. Пример оценки СПМ методами Уэлча и Берга для зашумленного двухчастотного сигнала

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

MATLAB-функция построения спектрограмм — specgram

Для визуализации короткооконного БПФ служит функция построения спектрограмм specgram. Это очень мощное и современное средство визуализации спектра, аппаратно реализованное в анализаторах спектра реального времени фирмы Tektronix [2, 3]. Спектрограмма представляется зависимостью амплитуды спектральных составляющих БПФ, вычисляемого в перемещающемся окне, от момента времени, задающего положение окна. Фактически спектрограмма строится в плоскости частота-время, а амплитуда каждой спектральной составляющей определяет цвет построения каждой точки спектрограммы. При построении спектрограммы используется функциональная окраска, иногда с применением шкалы цветов.

Данная функция имеет ряд форм записи. В простейшем случае B = specgram (a) вычисляется спектрограмма сигнала, с отчетами в векторе ж. При этом ряд параметров используется по умолчанию:

• 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.

Рис. 30. Спектрограмма звуковых колебаний

Алгоритм вычисления спектрограмм содержит три характерных шага:

1. Разбивка ж на перекрывающиеся блоки, на каждый из которых накладывается окно.
2. Выполнение nttf-точечного БПФ для соответствующего отрезка времени, что создает соответствующий столбец матрицы B, после чего окно перемещается на число точек, равное (length (window) — numoverlap). Если число точек БПФ превышает количество отсчетов в окне, то перед выполнением БПФ блок дополняется нулями.
3. При вещественных компонентах ж спектрограмма строится для положительных частот, и матрица B содержит при nfft четном nfft/2+1 строк, а при nfft нечетном — (nfft+1)/2 строк и k = fix((n-numoverlap)/ /(length(window)-numoverlap)) столбцов.

Подобные спектрограммы широко применяются в электроакустике, поскольку позволяют создавать красочный графический образ звуковых колебаний, в котором опытный взгляд может подметить множество тонких особенностей анализируемых звуков. Как шутку, в которой много правды, приведем в заключение спектрограмму колебания из файла vcosig (рис. 31), входящего в набор демонстрационных файлов системы МАTLАВ. Эта спектрограмма построена с помощью команд:

load vcosig; specgram(vcosig,[],Fs)

Рис. 31. Спектрограмма сигнала vcosig

Обработка в МАTLАВ спектрометрических данных в пакете биоинформатики

Одни из самых сложных спектров — спектры, получаемые масс-спектрометрами при исследовании состава веществ и биологических объектов. Их обработка является одним из основных методов контроля состава веществ в спектроскопии и спектрометрии, а порою и обнаружения новых веществ в микроскопических количествах. Ныне это крайне актуально при борьбе с терроризмом, когда нужно обнаруживать и распознавать малые дозы взрывчатых веществ и вещественные доказательства преступлений. Средства для работы с микромассивами веществ и обработки масс-спектрометрической информации сосредоточены в новом пакете расширения системы 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/Z)');ylabel('Относительная интенсивность')
title('Загрузка и построение спектров масс-спектрометрв')

Рис. 32. Построение четырех спектров масс-спектрометра

Спектры легко различаются по цвету. По пикам их кривых можно судить о составе веществ, для которых они получены. Для этого обычно требуется обработка спектров, например, перевыборка отчетов, коррекция базовой линии, поиск и сортировка пиков спектральных линий и т. д. Все эти возможности реализованы функциями пакета расширения Bioinformatics Toolbox. Следующий пример строит спектр с базовой линией (красной на рис. 33):

YB = msbackadj(MZ,Y,'WINDOWSIZE',500,'QUANTILE',0.20,
'SHOWPLOT',2);

Рис. 33. Построение спектра масс-спектрометра с базовой линией

Калибровка масс-спектрометра ведет к изменениям отношения между наблюдаемым M/Z вектором и истинным временем пролета ионов. Погрешности могут накапливаться при повторных экспериментах. Для устранения этого можно использовать функцию msalign, что позволит стандартизировать значения M/Z. В следующем примере задан массив спектрограммы и весов, используемых для осуществления алгоритма сглаживания:

PP = [3991.4 4598 7964 9160]; W = [60 100 60 100];

Теперь построим карту спектра (по сути, спектрограмму) до выравнивания, используя команды:

msheatmap(MZ,YB,'MARKERS',P,'LIMIT',[3000 10000])
title('После выравнивания (Alignment)')

Следующие команды обеспечивают построение карты спектра после процедуры выравнивания (рис. 34):

msheatmap(MZ,YB,'MARKERS',P,'LIMIT',[3000 10000])
title('Before Alignment')

Рис. 34. Карта спектра (спектрограмма) после выравнивания

Для уменьшения погрешностей при обработке спектра используется операция его нормализации. Интенсивность каждого сигнала приводится к некоторому значению, например 100. При этом можно исключить некоторые области спектра. Для нормализации используется функция msnorm.

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

Для выявления пиков спектрограмм можно использовать функцию численного дифференцирования МАTLАВ 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('Относительная интенсивность');
title('Автоматическая селекция пиков спектра')

Рис. 35. Спектр после обработки

Пакет имеет основанный на графическом интерфейсе пользователя GUI вьювер спектров масс-спектрометров. Вот пример его вывода функцией msviewer:

msviewer (MZ,YS,'MARKERS',MZ(h),'GROUP',1:4)

Спектры во вьювер могут загружаться из рабочего пространства MATLAB или из файлов формата CVS.

В последнее время помимо Фурье-спектрограмм сигналов стали использоваться вейв-лет-спектрограммы, которые имеют ряд преимуществ при анализе тонких особенностей сигналов. Техника вейвлет-преобразования сигналов и ее инструментальные средства (включая вейвлет-спектрограммы) описана в [5, 10]. Там же детально описан пакет расширения системы MATLAB по вейвле-там — Wavelet Toolbox. Применение техники вейвлет-преобразований (и Фурье-преобразований) при анализе реальных осциллограмм от современных цифровых осциллографов описано в [7].