Применение преобразования Фурье для анализа вибрационных сигналов.

Андрей Шелковенко. Один из разработчиков и основатель компании Vibromera.
Перевод статьи может содержать неточности.

1. Преобразование Фурье и спектр сигнала
   Во многих случаях задача получения (вычисления) спектра сигнала выглядит следующим образом. Имеется АЦП, который с частотой дискретизации Fd преобразует непрерывный сигнал, поступающий на его вход в течение времени T, в цифровые выборки - N штук. Затем этот массив выборок подается в некоторую программу (например FourierScope), который выводит N/2 некоторых числовых значений.
   Чтобы проверить правильность работы программы, формируем массив образцов как сумму двух sin(10*2*pi*x)+0,5*sin(5*2*pi*x) и подаем его в программу. Программа нарисовала следующее:

 

 

На графике спектра присутствуют две гармоники - 5 Гц с амплитудой 0,5 В и 10 Гц с амплитудой 1 В, все как в формуле исходного сигнала. Все в порядке, грограмма работает корректно.

Это означает, что если мы подадим на вход АЦП реальный сигнал, состоящий из смеси двух синусоид, то получим аналогичный спектр, состоящий из двух гармоник.

Итак, наш настоящий измеренный сигнал продолжительностью 5 сек., оцифрованные АЦП, т.е. представленные по отдельности образцы, имеет дискретный непериодический спектр.
С математической точки зрения - сколько ошибок в этой фразе?

Теперь попробуем измерить тот же сигнал в течение 0,5 с.

 

 

Что-то здесь не так! Гармоника на 10 Гц рисуется нормально, а вместо гармоники на 5 Гц появляются какие-то непонятные гармоники.

В Интернете пишут, что нужно добавить нули в конец образца и спектр будет рисоваться нормально.

 

 

Дело совсем не в этом. Мне придется разобраться с теорией. Давайте перейдем к википедия - источник знаний.

2. Непрерывная функция и ее представление в виде ряда Фурье
Математически наш сигнал длительностью T секунд - это некоторая функция f(x), заданная на интервале {0, T} (в данном случае X - это время). Такая функция всегда может быть представлена в виде суммы гармонических функций (синуса или косинуса) вида:

k - номер тригонометрической функции (номер гармонической составляющей, номер гармоники)
T - отрезок, на котором определяется функция (длительность сигнала)
Ak - амплитуда k-й гармонической составляющей,
θk - начальная фаза k-й гармонической составляющей
Что значит "представить функцию в виде суммы ряда"? Это значит, что, сложив значения гармонических составляющих ряда Фурье в каждой точке, мы получим значение нашей функции в этой точке.
(Более строго говоря, среднее квадратическое отклонение ряда от функции f(x) будет стремиться к нулю, но, несмотря на сходимость среднего квадрата, ряд Фурье функции, вообще говоря, не обязательно сходится к ней точечно. )
Этот ряд также можно записать в виде:

 

 

 

где k-я комплексная амплитуда.

 

или

 

 

 

Связь между коэффициентами (1) и (3) выражается следующими формулами:

 

 

 

 

Обратите внимание, что все эти три представления ряда Фурье совершенно эквивалентны. Иногда при работе с рядами Фурье удобнее вместо синусов и косинусов использовать экспоненты мнимого аргумента, то есть использовать преобразование Фурье в комплексной форме. Но нам удобнее пользоваться формулой (1), где ряд Фурье представлен в виде суммы косинусов с соответствующими амплитудами и фазами. В любом случае, неверно говорить, что результатом преобразования Фурье реального сигнала будут комплексные гармонические амплитуды. Как правильно сказано в Вики, "преобразование Фурье (ℱ) - это операция, которая отображает одну функцию вещественной переменной на другую функцию также вещественной переменной".

 

Итог:
Математической основой спектрального анализа сигналов является преобразование Фурье.

Преобразование Фурье позволяет представить непрерывную функцию f(x) (сигнал), заданную на интервале {0, T}, в виде суммы бесконечного числа (бесконечного ряда) тригонометрических функций (синуса и/или косинуса) с определенными амплитудами и фазами, также рассматриваемых на интервале {0, T}. Такой ряд называется рядом Фурье.

Отметим еще несколько моментов, понимание которых необходимо для правильного применения преобразования Фурье к анализу сигналов. Если мы рассмотрим ряд Фурье (сумму синусоид) по всей оси X, то увидим, что за пределами интервала {0, T} функция ряда Фурье будет периодически повторять нашу функцию.

Например, на графике рис. 7 исходная функция задана на интервале {-T\2, +T\2}, а ряд Фурье представляет собой периодическую функцию, заданную на всей оси x.

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

Таким образом:

Наша исходная функция - это непрерывная непериодическая функция, определенная на некотором отрезке длины T.
Спектр этой функции дискретен, то есть представлен в виде бесконечного ряда гармонических составляющих - ряда Фурье.
В самом деле, ряд Фурье определяет некоторую периодическую функцию, которая совпадает с нашей функцией на интервале {0, T}, но для нас эта периодичность не существенна.

Следующий.

Периоды гармонических составляющих кратны интервалу {0, T}, на котором определена исходная функция f(x). Другими словами, периоды гармоник кратны длительности измерения сигнала. Например, период первой гармоники в ряду Фурье равен интервалу T, на котором определена функция f(x). Период второй гармоники в ряду Фурье равен интервалу T/2. И так далее (см. рисунок 8).

Соответственно, частоты гармонических составляющих кратны 1/T. То есть, частоты гармонических составляющих Fk равны Fk= k\T, где k имеет значения от 0 до ∞, например, k=0 F0=0; k=1 F1=1\T; k=2 F2=2\T; k=3 F3=3\T;..... Fk= k\T (при нулевой частоте постоянная составляющая).

Пусть наша исходная функция - сигнал, записанный за время T=1 сек. Тогда период первой гармоники будет равен длительности нашего сигнала T1=T=1 сек, а частота гармоники равна 1 Гц. Период второй гармоники будет равен длительности нашего сигнала, деленной на 2 (T2=T/2=0,5 сек.), а частота равна 2 Гц. Для третьей гармоники T3=T/3 сек. и частота равна 3 Гц. И так далее.

Шаг между гармониками в этом случае составляет 1 Гц.

Так, сигнал длительностью 1 с может быть разложен на гармонические составляющие (для получения спектра) с частотным разрешением 1 Гц.
Для увеличения разрешения в 2 раза до 0,5 Гц необходимо увеличить длительность измерения в 2 раза до 2 с. 10-секундный сигнал может быть разложен на гармонические составляющие (спектр) с частотным разрешением 0,1 Гц. Других способов увеличить частотное разрешение не существует.

Существует способ искусственного увеличения длительности сигнала путем добавления нулей в массив выборок. Но это не увеличивает реальное частотное разрешение.

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

Обычная схема измерения и оцифровки сигнала выглядит следующим образом.

Измерительный преобразователь -- Нормализатор сигнала -- АЦП -- Компьютер
(Рис.9 Схема измерительного канала)

Сигнал с измерительного преобразователя поступает на АЦП в течение времени Т. Показания сигнала (выборка), полученные за время Т, передаются на компьютер и сохраняются в памяти.

Какие требования предъявляются к параметрам оцифровки сигнала? Устройство, преобразующее входной аналоговый сигнал в дискретный код (цифровой сигнал), называется аналого-цифровым преобразователем (АЦП) (© Wiki).

Одним из основных параметров АЦП является максимальная частота дискретизации - частота дискретизации непрерывного во времени сигнала. Частота дискретизации измеряется в герцах. ((© Wiki))

Согласно теореме Котельникова, если непрерывный сигнал имеет спектр, ограниченный частотой Fmax, то он может быть полностью и однозначно восстановлен из его дискретных выборок, взятых через интервалы времени T = 1/2*Fmax, то есть с частотой Fd ≥ 2*Fmax, где Fd - частота дискретизации; Fmax - максимальная частота спектра сигнала. Другими словами, частота оцифровки сигнала (частота дискретизации АЦП) должна быть как минимум в 2 раза выше максимальной частоты сигнала, который мы хотим измерить.

А что произойдет, если мы будем брать образцы с меньшей частотой, чем того требует теорема Котельникова?

В этом случае возникает эффект "алиасинга" (он же стробоскопический эффект, эффект муара), при котором высокочастотный сигнал после оцифровки превращается в низкочастотный, которого на самом деле не существует. На рис. 11 красная синусоида высокой частоты - это реальный сигнал. Синяя синусоида более низкой частоты - фиктивный сигнал, возникающий из-за того, что за время дискретизации успевает пройти более половины периода высокочастотного сигнала.

 

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

Для вычисления спектра сигнала по его дискретным выборкам используется дискретное преобразование Фурье (ДПФ). Заметим еще раз, что спектр дискретного сигнала "по определению" ограничен частотой Fmax, меньшей половины частоты дискретизации Fd. Поэтому спектр дискретного сигнала можно представить в виде суммы конечный количество гармоник, в отличие от бесконечной суммы для ряда Фурье непрерывного сигнала, спектр которого может быть неограниченным. Согласно теореме Котельникова, максимальная частота гармоники должна быть такой, чтобы на нее приходилось не менее двух отсчетов, поэтому число гармоник равно половине числа отсчетов дискретного сигнала. То есть, если в выборке N выборок, то количество гармоник в спектре будет N/2.

Теперь рассмотрим дискретное преобразование Фурье (ДПФ).

Сравнение с рядом Фурье

 

Как видим, они совпадают, за исключением того, что время в БПФ дискретно и количество гармоник ограничено N/2, что в два раза меньше числа отсчетов.

Формулы ДПФ записываются в безразмерных целочисленных переменных k, s, где k - количество выборок сигнала, s - количество спектральных компонент.
Величина s показывает количество полных гармонических колебаний за период T (длительность измерения сигнала). Дискретное преобразование Фурье используется для численного, т.е. "компьютерного", нахождения амплитуд и фаз гармоник.

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

 

 

Как видно из рис. 12, функция f(x) является периодической с периодом T0. Однако из-за того, что длина измерительного образца T не равна периоду функции T0, функция, полученная в виде ряда Фурье, имеет разрыв в точке T. В результате спектр этой функции будет содержать большое количество высокочастотных гармоник. Если бы длительность измерения образца T совпадала с периодом функции T0, то спектр, полученный после преобразования Фурье, содержал бы только первую гармонику (синусоиду с периодом, равным длительности образца), поскольку функция f(x) является синусоидой.

Другими словами, программа DFT "не знает", что наш сигнал - это "кусочек синусоиды", но пытается представить в виде ряда периодическую функцию, которая имеет прерывистость, обусловленную прерывистостью отдельных кусочков синусоиды.

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

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

 

 

При меньшей длительности изображение будет выглядеть "хуже":

 

 

 

 

На практике бывает трудно понять, где "реальные компоненты", а где "артефакты", вызванные несоответствием периодов компонентов и длительностей дискретизации сигнала или "скачками и разрывами" в форме сигнала. Конечно, слова "реальные компоненты" и "артефакты" взяты в кавычки не случайно. Наличие множества гармоник на графике спектра не означает, что наш сигнал действительно состоит из них. Это все равно что считать, что число 7 "состоит" из чисел 3 и 4. Число 7 можно представить как сумму 3 и 4 - это верно.

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

Некоторые результаты
1. Реальный измеряемый сигнал длительностью T сек., оцифрованный АЦП, т.е. представленный набором дискретных выборок (N штук), имеет дискретный непериодический спектр, представленный набором гармоник (N/2 штук).

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

3. Сигнал, измеренный в момент времени T, определяется только в момент времени T. Что было до того, как мы начали измерять сигнал, и что будет после - науке неизвестно. И в нашем случае это неинтересно. БПФ ограниченного по времени сигнала дает его "реальный" спектр, в том смысле, что при определенных условиях позволяет вычислить амплитуду и частоту его составляющих.