Электронная библиотека
по теме квалификационной работы магистра
Кикалова Владимира Сергеевича
"Разработка экспертной системы для распознования причин повреждений элементов
металлургических машин"
Руководитель: доц. каф. МОЗЧМ, к.т.н. Сидоров Владимир Анатольевич
ПРЕОБРАЗОВАНИЕ
ФУРЬЕ И КЛАССИЧЕСКИЙ ЦИФРОВОЙ СПЕКТРАЛЬНЫЙ АНАЛИЗ. Введение Спектральный анализ - один из методов обработки сигналов, который позволяет охарактеризовать частотный состав измеряемого сигнала. Преобразование Фурье является математической основой, которая связывает временной или пространственный сигнал (или же некоторую модель этого сигнала) с его представлением в частотной области. Важную роль в спектральном анализе играют методы статистики, поскольку сигналы, как правило, имеют случайный характер или зашумлены при распространении или измерении. Если бы основные статистические характеристики сигнала были точно известны, или их можно было определить по конечному интервалу этого сигнала, то спектральный анализ представлял бы собой отрасль "точной науки". Однако, в действительности по отрезку сигнала можно получить только оценку его спектра. Поэтому практика спектрального анализа - некое ремесло (или искусство?) достаточно субъективного характера. Различие между спектральными оценками, получаемыми в результате обработки одного и того же отрезка сигнала разными методами, можно объяснить различием допущений, принятых относительно данных, различными способами усреднения и т .п. Если априори характеристики сигнала не известны, нельзя сказать какие из оценок лучше. Преобразование
Фурье - математическая основа спектрального анализа ,
(1) которое
идентифицирует частоты и амплитуды тех комплексных синусоид (экспонент),
на которые разлагается некоторое произвольное колебание.
Менее ограничительное достаточное условие - конечность энергии сигнала
а функция sinc - выражением
Функция
отсчетов во временной области определяется выражением (7)
Таблица 1. Основные свойства НВПФ и функции
Еще одно важное свойство устанавливается теоремой Парсеваля для двух функций g(t) и h(t):
Если
положить g(t) = h(t), то теорема Парсеваля сводится к теореме для энергии . (9) Выражение (9) - это, в сущности, просто формулировка закона сохранения энергии в двух областях (временной и частотной). В (9) слева стоит полная энергия сигнала, таким образом, функция
описывает распределение энергии по частоте для детерминированного сигнала h(t) и поэтому называется спектральной плотностью энергии (СПЭ). С помощью выражений
можно вычислить амплитудный и фазовый спектры сигнала h( t ). Операции дискретизации и взвешивания В следующем разделе мы введем дискретно-временной ряд Фурье (ДВРФ) или иначе дискретное преобразование Фурье (ДПФ) как частный случай непрерывно-временного преобразования Фурье (НВПФ) с использованием двух базовых операций обработки сигналов - взятия отсчетов (дискретизации) и взвешивания с помощью окна. Здесь рассмотрим влияние этих операций на сигнал и его преобразование. В таблице 2 перечислены функции, с помощью которых осуществляется взвешивание и дискретизация . При равномерных отсчетах с интервалом T секунд частота отсчетов F равна 1 /T Гц. Заметим, что взвешивающая функция и функция отсчетов во временной области обозначаются соответственно TW (time windowing) и TS (time sampling), а в частотной области - FW (frequency windowing) и FS (frequency sampling).
Предположим,
что берутся отсчеты непрерывного действительного сигнала x(t) c ограниченным
спектром, верхняя частота которого равна F0. НВПФ действительного сигнала
- это всегда симметричная функция с полной шириной 2F0, см. рис.1.
Рис.1
- иллюстрация теоремы отсчетов во временной области для действительного
сигнала с ограниченным спектром:
Свертка
X(f) c преобразованием Фурье функции отсчетов F {TS}=Y1/T( f ) просто
периодически продолжает X(f) с частотным интервалом 1/T Гц. Поэтому XS(f)
представляет собой периодически продолженный спектр X(f). В общем случае
отсчеты в одной области (например, временной) приводят к периодическому
продолжению в области преобразования (например, частотной). Если частота
отсчетов выбрана достаточно низкой (F < 2Fo ), то периодически продолженные
спектры будут перекрываться с соседними. Это перекрытие носит название
эффекта наложения в частотной области.
В результате (см. Рис. 1 д ) восстанавливается исходное преобразование Фурье. Используя теоремы о свертке во временной и частотной областях, получаем . (15) Выражение
(15) представляет собой математическую запись теоремы отсчетов во временной
области ( теоремы Уиттекера, Котельникова, Шеннона - УКШ ), которая
утверждает, что с помощью интерполяционной формулы (15) действительный
сигнал с ограниченным спектром может быть точно восстановлен по бесконечному
числу известных временных отсчетов, взятых с частотой F і 2F0. Дуальной
к теореме (15) является теорема отсчетов в частотной области для
сигналов с ограниченной длительностью. , (16) а соответствующие преобразования - выражениями
Таким образом, НВПФ X(f) некоторого сигнала с ограниченной длительностью может быть однозначно восстановлено по эквидистантным отсчетам спектра такого сигнала, если выбранный интервал отсчетов по частоте удовлетворяет условию F1/2T0 Гц, где T0 - длительность сигнала. Соотношения между непрерывными и дискретными преобразованиями Пара преобразований для обычного определения дискретного преобразования Фурье (ДПФ) N-точечной временной последовательности x[n] и соответствующей ей N-точечной последовательности преобразования Фурье X[k] дается выражениями ,
(18) Чтобы по отсчетам данных получить спектральные оценки в соответствующих единицах измерения энергии или мощности, запишем дискретно-временной ряд Фурье (ДВРФ), который можно рассматривать как некоторую аппроксимацию непрерывно-временного преобразования Фурье (НВПФ), основанную на использовании конечного числа отсчетов данных:
Для того, чтобы показать характер соответствия ДВРФ (дискретные функции и во временной и в частотной областях) и НВПФ (непрерывные функции во временной и в частотной областях), нам потребуется последовательность из четырех линейных коммутативных операций: взвешивания во временной и частотной областях и взятия отсчетов или дискретизации как во временной, так и в частотной областях. Если операция взвешивания выполняется в одной из этих областей, то, согласно теореме свертки, ей будет соответствовать выполнение операции фильтрации (свертки) в другой области с функцией sinc. Точно также, если в одной области выполняется дискретизация, то в другой выполняется операция периодического продолжения. Так как взвешивание и взятие отсчетов являются линейными и коммутативными операциями, то возможны различные способы их упорядочения, дающие одинаковый конечный результат при различных промежуточных результатах. На рис.2 показаны две возможные последовательности выполнения этих четырех операций. Рис.
2. Две возможные последовательности из двух операций взвешивания и двух
операций взятия отсчетов, связывающие НВПФ и ДВРФ: FW - применение окна
в частотной области; TW - применение окна во временной области; FS - взятие
отсчетов в частотной области; TS - взятие отсчетов во временной области.
(22) (23)
( 24 ) Таким
образом, ДВПФ - это просто Z-преобразование, вычисленное на единичной
окружности и умноженное на T.
где k=-N/2, . . . ,N/2-1
где
n=0, . . . ,N-1 , ,
(29) и характеризует энергию последовательности из N отсчетов данных. Обе последовательности x[n] и X[k] периодичны по модулю N, поэтому (28) можно записать в форме , (30) где 0 n N. Множитель T в (27) - (30) необходим для того, чтобы (27) и (28) являлись в действительности аппроксимацией интегрального преобразования в области интегрирования .(31) Дополнение нулями С помощью процесса, называемого дополнением нулями, дискретно-временной ряд Фурье может быть изменен для интерполяции между N значениями исходного преобразования. Пусть имеющиеся отсчеты данных x[0],...,x[N-1] дополнены нулевыми значениями x[N],...X[2N-1]. ДВРФ этой дополненной нулями 2N-точечной последовательности данных будет определяться выражением
(32) где
верхний предел суммы справа изменен в соответствии с наличием нулевых
данных. Пусть k=2m, так что ,
(33) где m=0,1,...,N-1, определяет четные значения X[k]. Отсюда видно, что при четных значениях индекса k 2N-точечный дискретно-временной ряд Фурье сводится к N-точечному дискретно-временному ряду. Нечетные значения индекса k соответствуют интерполированным значениям ДВРФ, расположенным между значениями исходного N-точечного ДВРФ. По мере того, как все большее число нулей добавляется в исходную N-то-чечную последовательность, можно получить еще большее число интерполированных данных. В предельном случае бесконечного числа вводимых нулей ДВРФ может рассматриваться как дискретно-временное преобразование Фурье N-то-чечной последовательности данных:
Преобразование
(34) соответствует узлу 6 на рис.2.
Рис.3.
Интерполяция за счет дополнения нулями:
Можно показать, что эквивалентная длительность окна (или любого другого сигнала) и эквивалентная ширина полосы его преобразования являются взаимно обратными величинами: TeBe=1. Быстрое преобразование Фурье Быстрое
преобразование Фурье (БПФ) - это не еще одна разновидность преобразования
Фурье, а название целого ряда эффективных алгоритмов, предназначенных
для быстрого вычисления дискретно-временного ряда Фурье. Основная
проблема, возникающая при практической реализации ДВРФ, заключена в большом
количестве вычислительных операций, пропорциональном N2. Хотя еще задолго
до появления компьютеров было предложено несколько эффективных вычислительных
схем, позволяющих существенно сократить число вычислительных операций,
настоящую революцию произвела публикация в 1965 году статьи Кули (Cooly)
и Тьюки (Tukey) c практическим алгоритмом быстрого (число операций Nlog2N)
вычисления ДВРФ. После этого было разработано множество вариантов, усовершенствований
и дополнений основной идеи, составивших класс алгоритмов, известных под
названием быстрого преобразования Фурье. Основная идея БПФ - деление N-точечного
ДВРФ на два и более ДВРФ меньшей длины, каждый из которых можно вычислить
отдельно, а затем линейно просуммировать с остальными, с тем чтобы получить
ДВРФ исходной N-точечной последовательности.
, (35) где величина WN=exp(-j2/N) носит название поворачивающего множителя (здесь и далее в этом разделе период выборки T=1). Выделим из последовательности x[n] элементы с четными и нечетными номерами
Но
так как то , (37) где каждое из слагаемых является преобразованием длины N/2 (38) Заметим, что последовательность (WN/2)nk периодична по k с периодом N/2. Поэтому, хотя номер k в выражении (37) принимает значения от 0 до N-1, каждая из сумм вычисляется для значений k от 0 до N/2-1. Можно оценить число комплексных операций умножения и сложения, необходимых для вычисления преобразования Фурье в соответствии с алгоритмом (37)-(38). Два N/2-точечных преобразования Фурье по формулам (38) предполагают выполнение 2(N/2)2 умножений и приблизительно столько же сложений. Объединение двух N/2-точечных преобразований по формуле (37) требует еще N умножений и N сложений. Следовательно, для вычисления преобразования Фурье для всех N значений k необходимо произвести по N+N2/2 умножений и сложений. В то же время прямое вычисление по формуле (35) требует по N2 умножений и сложений. Уже при N>2 выполняется неравенство N+N2/2 < N2 , и, таким образом, вычисления по алгоритму (37)-(38) требуют меньшего числа математических операций по сравнению с прямым вычислением преобразования Фурье по формуле (35). Так как вычисление N-точечного преобразования Фурье через два N/2-точечных приводит к экономии вычислительных операций, то каждое из N/2-точечных ДПФ следует вычислять путем сведения их к N/4-точечным преобразованиям: ,
(39)
(41) где "одноточечные преобразования" X1[k,p] представляют собой просто отсчеты сигнала x[n]: X1[k,q]
= x[q]/N, q=0,1,...,N-1.
(42) В итоге можно записать алгоритм БПФ, получивший по понятным причинам название алгоритма с прореживанием по времени : X2[k,p] = (x[p] + Wk2x[p+N/2]) / N, где k=0,1, p=0,1,...,N/2 -1; X2N/M[k,p] =XN/M[k,p] + Wk2N/MXN/M[k,p+M/2], где k=0,1,...,2N/M -1, p=0,1,...,M/2 -1; ..X[k] = XN[k] =XN/2[k,0] + WkNXN/2[k,1], (43) где k=0,1,...,N-1 На
каждом этапе вычислений производится по N комплексных умножений и сложений.
А так как число разложений исходной последовательности на подпоследовательности
половинной длины равно log2N, то полное число операций умножения-сложения
в алгоритме БПФ равно Nlog2N. При больших N имеет место существенная
экономия вычислительных операций по сравнению с прямым вычислением ДПФ.
Например, при N = 210 = 1024 число операций уменьшается в 117
раз. Случайные процессы и спектральная плотность мощности Дискретный случайный процесс x[n,i] можно рассматривать как некоторую совокупность, или ансамбль, действительных или комплексных дискретных временных (или пространственных) последовательностей, каждую из которых можно было бы наблюдать как результат проведения некоторого эксперимента (n - временной индекс, i - номер наблюдения). Последовательность, полученную в результате одного из наблюдений будем обозначать x[n]. Операцию усреднения по ансамблю (т.е. статистического усреднения) будем обозначать посредством оператора <>. Таким образом, <x[n]> - среднее значение случайного процесса x[n] в момент времени n. Автокорреляция случайного процесса в два различных момента времени n1 и n2 определяется выражением rxx[n1,n2]=<x[n1]x*[n2]>. Случайный
процесс называется стационарным в широком смысле, если его среднее
значение постоянно (не зависит от времени), а автокорреляция зависит только
от разности индексов времени m=n1-n2 (временного сдвига или задержки между
отсчетами). Таким образом, стационарный в широком смысле дискретный случайный
процесс x[n] характеризуется постоянным средним значением <x>=<x[n]>
и автокорреляционной последовательностью (АКП) rxx[m] = < x[n+m]x*[n] >. (44) Отметим следующие свойства АКП: rxx[0] |rxx[m]| , rxx[-m] = r*xx[m] , (45) которые
справедливы при всех m. . (46) СПМ, ширина которой полагается ограниченной значениями ±1/2T Гц, является периодической функцией частоты с периодом 1/T Гц. Функция СПМ описывает распределение мощности случайного процесса по частоте. Для подтверждения избранного для нее названия рассмотрим обратное ДВПФ (47) вычисляемое при m=0 (48) Автокорреляция при нулевом сдвиге характеризует среднюю мощность случайного процесса. Согласно (48), площадь под кривой Pxx(f) характеризует среднюю мощность, поэтому Pxx(f) представляет собой функцию плотности (мощность на единицу измерения частоты), которая характеризует распределение мощности по частоте. Пару преобразований (46) и (47) часто называют теоремой Винера-Хинчина для случая дискретного времени. Поскольку rxx[-m]=r*xx[m], то СПМ должна быть строго действительной положительной функцией. Если АКП - строго действительная функция, то rxx[-m]=rxx[m] и СПМ можно записать в форме косинус-преобразования Фурье , что
означает также, что Pxx(f) = Pxx(-f), т.е. СПМ -
четная функция. . (49) Этот предел, если он существует, сходится к истинному среднему значению тогда и только тогда, когда дисперсия среднего по времени стремится к нулю, что означает выполнение следующего условия: . (50)
(51) Допущение об эргодичности позволяет не только ввести через усреднение по времени определения для среднего значения и автокорреляции, но также дать подобное определение и для спектральной плотности мощности . (52) Эта эквивалентная форма СПМ получается посредством статистического усреднения модуля ДВПФ взвешенной совокупности данных, поделенного на длину записи данных, для случая, когда число отсчетов увеличивается до бесконечности. Статистическое усреднение необходимо здесь потому, что ДВПФ само является случайной величиной, изменяющейся для каждой реализации x[n]. Для того, чтобы показать, что (52) эквивалентно теореме Винера-Хинчина, представим квадрат модуля ДВПФ в виде произведения двух рядов и изменим порядок операций суммирования и статистического усреднения: (53) Используя известное выражение , (54)
( 55 ) Заметим, что на последнем этапе вывода (55) использовалось допущение о том, что автокорреляционная последовательность "затухает", так что . (56) Взаимосвязь
двух определений СПМ (46) и (52) наглядно показывает диаграмма, представленная
на рисунке 4. , ( 57) которая называется выборочным спектром.
Рис. 4. Взаимосвязь двух способов оценивания спектральной плотности мощности Периодограммный метод спектрального оценивания Выше мы ввели два формальных эквивалентных метода определения спектральной плотности мощности (СПМ). Косвенный метод основан на использовании бесконечной последовательности данных для расчета автокорреляционной последовательности, преобразование Фурье которой дает искомую СПМ. Прямой метод определения СПМ основан на вычислении квадрата модуля преобразования Фурье для бесконечной последовательности данных с использованием соответствующего статистического усреднения. СПМ, полученная без такого усреднения оказывается неудовлетворительной, поскольку среднеквадратичная ошибка такой оценки сравнима с ее средним значением. Сейчас мы рассмотрим методы усреднения, обеспечивающие получение гладких и статистически устойчивых спектральных оценок по конечному числу отсчетов. Оценки СПМ, основанные на прямом преобразовании данных и последующем усреднении, получили название периодограмм. Оценки СПМ, для получения которых по исходным данным сначала формируются корреляционные оценки, получили название коррелограммных. При использовании любого метода оценивания СПМ пользователю приходится принимать множество компромиссных решений, с тем чтобы по конечному количеству отсчетов получать статистически устойчивые спектральные оценки с максимально возможным разрешением. К этим компромиссам относятся, в частности, выбор окна для взвешивания данных и корреляционных оценок и таких параметров усреднения во временной и в частотной областях, которые позволяют сбалансировать требования к снижению уровня боковых лепестков, обусловленных взвешиванием, выполнению эффективного усреднения и к обеспечению приемлемого спектрального разрешения. На рис. 5 приведена диаграмма, отображающая основные этапы периодограммного метода
Применение метода начинается со сбора N отсчетов данных, которые берутся с интервалом T секунд на отсчет с последующим (по желанию) этапом устранения тренда. Для того, чтобы получить статистически устойчивую спектральную оценку, имеющиеся данные необходимо разбить на перекрывающиеся (по возможности) сегменты и в последующем усреднить выборочные спектры, полученные по каждому такому сегменту. Параметры этого усреднения изменяются посредством соответствующего выбора числа отсчетов на сегмент (NSAMP) и числа отсчетов, на которое необходимо сдвинуть начало следующего сегмента (NSHIFT), см. рис. 6. Количество сегментов выбирается в зависимости от требуемой степени гладкости (дисперсии) спектральной оценки и требуемого спектрального разрешения. При малом значении параметра NSAMP получается больше сегментов, по которым будет производиться усреднение, а следовательно будут получаться оценки с меньшей дисперсией, но также и меньшим частотным разрешением. Увеличение длины сегмента (параметра NSAMP) повышает разрешение, естественно за счет увеличения дисперсии оценки из-за меньшего числа усреднений. Стрелка возврата на рис.5 указывает на необходимость нескольких повторных проходов по данным при различных длинах и количествах сегментов, что позволяет получить больше информации об исследуемом процессе. Рис.6. Разбиение данных на сегменты для вычисления периодограммы Окна Один из важных вопросов, который является общим для всех классических методов спектрального оценивания, связан с взвешиванием данных. Обработка с помощью окна используется для управления эффектами, обусловленными наличием боковых лепестков в спектральных оценках. Заметим, что имеющуюся конечную запись данных удобно рассматривать как некоторую часть соответствующей бесконечной последовательности, видимую через применяемое окно. Так последовательность наблюдаемых данных x0[n] из N отсчетов математически можно записать как произведение бесконечной последовательности x[n] и функции прямоугольного окна x0[n]=x[n]·rect[n].
X0(f)=X(f)*DN(f)
, где Функция DN(f), называемая дискретной функцией sinc, или ядром Дирихле, представляет собой ДВПФ прямоугольной функции. Преобразование наблюдаемой конечной последовательности является искаженной версией преобразования бесконечной последовательности. Влияние прямоугольного окна на дискретно-временную синусоиду с частотой f0 иллюстрирует рис.7. Рис.7. Иллюстрация смещения дискретно-временного преобразования Фурье вследствие просачивания из-за взвешивания данных.: а,в - исходная и взвешенная последовательности; б, г - их преобразования Фурье. Из
рисунка видно, что острые спектральные пики ДВПФ бесконечной синусоидальной
последовательности расширились за счет свертки с преобразованием окна.
Таким образом, минимальная ширина спектральных пиков взвешенной окном
последовательности определяется шириной главного лепестка преобразования
этого окна и не зависит от данных. Боковые лепестки преобразования окна
будут изменять амплитуды соседних спектральных пиков (иногда это явление
называют просачиванием). Поскольку ДВПФ - периодическая функция, то наложение
боковых лепестков от соседних периодов может привести к дополнительному
смещению. Увеличение частоты отсчетов позволяет ослабить эффект наложения
боковых лепестков. Аналогичные искажения будут, естественно, наблюдаться
и в случае несинусоидальных сигналов. Просачивание приводит не только
к появлению амплитудных ошибок в спектрах дискретных сигналов, но может
также маскировать присутствие слабых сигналов. Можно предложить ряд других
функций окна, применение которых позволяет снизить уровень боковых лепестков
по сравнению с тем, который имеется при использовании прямоугольного окна.
Снижение уровня боковых лепестков будет уменьшать смещение спектральной
оценки, однако это дается ценой расширения главного лепестка спектра окна,
что, естественно, приводит к ухудшению разрешения. Следовательно и здесь
должен выбираться какой-то компромисс между шириной главного лепестка
и уровнем боковых лепестков. Для оценки качества окон используется несколько
параметров. Традиционным показателем является ширина полосы главного лепестка
на уровне половинной мощности. В качестве второго показателя используется
эквивалентная ширина полосы, введенная выше. Два показателя используются
и для оценки характеристик боковых лепестков. Первый - их максимальный
уровень, второй - скорость спадания, характеризующая быстроту уменьшения
боковых лепестков по мере удаления от главного лепестка. В таблице 3 приведены
определения некоторых общеупотребительных дискретно-временных функций
окна, а в таблице 4 - их характеристики.
Таблица 4. Некоторые характеристики типичных N-точечных дискретно-временных окон
Ранее мы дали определение спектральной плотности мощности как дискретно-временное преобразование Фурье автокорреляционной последовательности . (46) Коррелограммный метод оценивания СПМ - это просто подстановка в выражение (46) конечной последовательности значений оценки автокорреляции (коррелограммы) вместо бесконечной последовательности неизвестных истинных значений автокорреляции. Подробнее о коррелограммном методе спектрального оценивания можно прочитать в [1]. Л и т е р а т у р а 1.
Рабинер Л., Гоулд Б. Теория и применение цифровой обработки сигналов.
М.:Мир, 1978. 2.
Марпл-мл. С.Л. Цифровой спектральный анализ и его приложения: Пер. с англ.
-М.: Мир, 1990. 3.
Гольдберг Л.М., Матюшкин Б.Д., Поляк М.Н., Цифровая обработка сигналов.-
М.: Радио и связь, 1990. 4. Отнес Р., Эноксон Л. Прикладной анализ временных рядов.- М.: Мир, 1982.
|