Застосування перетворення Фур'є для аналізу сигналів вібрації
Андрій Шелковенко. Один з розробників та засновник компанії Vibromera.
Переклад статті може містити неточності.
Перетворення Фур'є та спектр сигналу
У багатьох випадках завдання отримання (обчислення) спектр сигналу виглядає наступним чином. Є АЦП, який за допомогою дискретизації частота Fd перетворює неперервний сигнал, що надходить на його вхід протягом часу T, у цифрові відліки – N штук. Потім цей масив відліків подається до деякої програми (наприклад Фур'є-скоп), який виводить N/2 деяких числових значень.
Для перевірки коректності роботи програми сформуємо масив відліків у вигляді суми двох sin(10*2*pi*x)+0.5*sin(5*2*pi*x) і передамо його в програму. Програма намалювала наступне:
Рис.1 Графік часової функції сигналу
Рис.2 Графік спектру сигналу
Існує два гармоніки на графіку спектра – 5 Hz з амплітудою 0.5 V і 10 Hz з амплітудою 1 V, усе як у формулі вихідного сигналу. Усе гаразд, програма працює правильно.
Це означає, що якщо ми подамо на вхід АЦП реальний сигнал із суміші двох синусоїд, то отримаємо подібний спектр, що складається з двох гармонік.
Отже, наш справжній виміряний сигнал тривалістю 5 секундоцифровані АЦП, тобто представлені дискретними зразки, має дискретний неперіодичний спектру.
З математичної точки зору – скільки помилок у цій фразі?
Тепер спробуймо виміряти той самий сигнал протягом 0.5 sec.
Рис.3 Графік функції sin(10*2*pi*x)+0.5*sin(5*2*pi*x) для періоду вимірювання 0.5 сек
Рис.4 Спектр функції
Щось тут не так! Гармоніка на 10 Гц намальована нормально, а замість гармоніки на 5 Гц є якісь незрозумілі гармоніки.
В інтернеті пишуть, що треба додати нулі в кінець зразка і спектр намалюється нормально.
Рис.5 Ми додали нулі до вибірки до 5 сек
Рис.6. Отриманий спектр.
Це зовсім не те. Доведеться розібратися з теорією. Перейдімо до вікіпедія – джерела знань.
Неперервна функція та її представлення у вигляді ряду Фур'є
Математично наш сигнал тривалістю T секунд є деякою функцією f(x), заданою на інтервалі {0, T} (X в даному випадку - час). Таку функцію завжди можна представити у вигляді суми гармонічних функцій (синуса або косинуса) виду:
(1), де:
k - номер тригонометричної функції (номер гармонічної складової, номер гармоніки)
T – відрізок, на якому визначена функція (тривалість сигналу)
Ak - амплітуда k-ї гармонічної складової,
θk - початкова фаза k-ї гармонічної складової
Що означає "представити функцію у вигляді суми ряду"? Це означає, що додаючи значення гармонічних складових ряду Фур'є в кожній точці, ми отримуємо значення нашої функції в цій точці.
(Більш строго, середнє квадратичне відхилення ряду від функції f(x) буде прагнути до нуля, але, незважаючи на збіжність середнього квадратичного, ряд Фур'є функції не обов'язково, взагалі кажучи, збігається до неї точка за точкою. )
Цю серію також можна написати у формі:
(2),
де 
k-та комплексна амплітуда.
або
(3)
Зв'язок між коефіцієнтами (1) і (3) виражається наступними формулами:
Зверніть увагу, що всі ці три представлення ряду Фур’є цілком еквівалентні. Іноді при роботі з рядами Фур’є зручніше використовувати показники з уявним аргументом замість синусів і косинусів, тобто застосовувати перетворення Фур’є в комплексній формі. Але для нас зручно використовувати формулу (1), де ряд Фур’є подано як суму косинусів із відповідними амплітудами та фазами. Строго кажучи, перетворення Фур’є дійсного сигналу дійсно дає комплексні коефіцієнти (форма (3)): кожен коефіцієнт несе і амплітуду, і фазу своєї гармоніки. Для дійсного сигналу ці комплексні коефіцієнти мають спряжену (ермітову) симетрію — від’ємна частотна половина просто віддзеркалює додатну частотну половину і не несе додаткової інформації. Саме тому з комплексних коефіцієнтів ми завжди можемо перейти до дійсних невід’ємних амплітуд Ak і фаз θk формули (1) — і саме цей амплітудний спектр відображають програми аналізу.
Підсумок:
Математичною основою спектрального аналізу сигналів є перетворення Фур'є.
Перетворення Фур'є дозволяє представити неперервну функцію f(x) (сигнал), визначену на проміжку {0, T}, у вигляді суми нескінченної кількості (нескінченного ряду) тригонометричних функцій (синуса та/або косинуса) з певними амплітудами і фазами, також розглянутих на проміжку {0, T}. Такий ряд називається рядом Фур'є.
Відзначимо ще кілька моментів, розуміння яких необхідне для коректного застосування перетворення Фур'є до аналізу сигналів. Якщо ми розглянемо ряд Фур'є (суму синусоїд) на всій осі Х, то побачимо, що за межами інтервалу {0, T} функція ряду Фур'є буде періодично повторювати нашу функцію.
Наприклад, на графіку на рис. 7 вихідна функція визначена на інтервалі {-T\2, +T\2}, а ряд Фур'є представляє періодичну функцію, визначену на всій осі х.
Це пояснюється тим, що самі синусоїди є періодичними функціями, тому їхня сума також буде періодичною функцією.
Рисунок 7 Представлення неперіодичної функції джерела за допомогою ряду Фур'є
Таким чином:
Наша вихідна функція є неперервною, неперіодичною функцією, визначеною на деякому відрізку довжини T.
Спектр цієї функції є дискретним, тобто представлений у вигляді нескінченного ряду гармонійних складових - ряду Фур'є.
Насправді, ряд Фур'є визначає деяку періодичну функцію, яка збігається з нашою функцією на проміжку {0, T}, але для нас ця періодичність не є суттєвою.
Наступний.
Періоди гармонічних складових кратні інтервалу {0, T}, на якому визначена початкова функція f(x). Іншими словами, періоди гармонік кратні тривалості вимірювання сигналу. Наприклад, період першої гармоніки в ряді Фур'є дорівнює інтервалу T, на якому визначена функція f(x). Період другої гармоніки в ряді Фур'є дорівнює інтервалу T/2. І так далі (див. Рисунок 8).
Рис. 8 Періоди (частоти) гармонійних складових ряду Фур'є (тут T=2π)
Відповідно, частоти гармонічних складових кратні 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 Гц.
Щоб підвищити роздільну здатність удвічі — до 0,5 Гц, необхідно збільшити тривалість вимірювання удвічі — до 2 секунд. 10-секундний сигнал можна розкласти на гармонічні складові (спектр) із частотною роздільною здатністю 0,1 Гц. Інших способів підвищити частотну роздільну здатність не існує. Ви можете дослідити цю залежність за допомогою нашого Калькулятор роздільної здатності FFT.
Існує спосіб штучного збільшення тривалості сигналу шляхом додавання нулів до масиву відліків. Але це не збільшує реальну частотну роздільну здатність.
Дискретні сигнали та дискретне перетворення Фур’є
З розвитком цифрових технологій змінилися способи зберігання вимірювальних даних (сигналів). Якщо раніше сигнал можна було записати на магнітофон і зберігати на стрічці в аналоговому вигляді, то зараз сигнали оцифровуються і зберігаються у файлах в пам'яті комп'ютера у вигляді набору чисел (відліків).
Звичайна схема вимірювання та оцифрування сигналу виглядає наступним чином.
Вимірювальний перетворювач —- Нормалізатор сигналу —- АЦП —– Комп'ютер
(Рис.9 Схема вимірювального каналу)
Сигнал від вимірювального перетворювача надходить до АЦП протягом періоду часу T. Отримані за час T відліки сигналу (дискретизація) передаються до комп'ютера і зберігаються в пам'яті.
Рис.10 Оцифрований сигнал – N відліків, отриманих за час T
Які вимоги висуваються до параметрів оцифрування сигналу? Пристрій, який перетворює вхідний аналоговий сигнал у дискретний код (цифровий сигнал), називається аналого-цифровим перетворювачем (АЦП) (© Wiki).
Одним з основних параметрів АЦП є максимальна частота дискретизації – частота дискретизації неперервного в часі сигналу. Частота дискретизації вимірюється в герцах. ((© Wiki))
According to Kotelnikov’s theorem, if a continuous signal has a spectrum limited by the frequency Fmax, it can be fully and uniquely reconstructed from its discrete samples taken at time intervals Δt ≤ 1/(2*Fmax), ie with a sampling frequency Fd ≥ 2*Fmax, where Fd – sampling frequency; Fmax – the maximum frequency of the signal spectrum. In other words, the frequency of signal digitization (sampling frequency of ADC) must be at least twice the maximum frequency of the signal we want to measure.
А що станеться, якщо брати відліки з меншою частотою, ніж вимагає теорема Котельникова?
У цьому випадку виникає “аліасингуУ цьому випадку виникає ефект "аліасингу" (він же стробоскопічний ефект, ефект муару), при якому високочастотний сигнал після оцифрування перетворюється на низькочастотний сигнал, якого насправді не існує. На рис. 11 червона синусоїда високої частоти є реальним сигналом. Синя синусоїда нижчої частоти - це фіктивний сигнал, який виникає через те, що за час дискретизації встигає пройти більше половини періоду високочастотного сигналу.
Рис. 11. Поява паразитного низькочастотного сигналу при недостатньо високій частоті дискретизації
Щоб уникнути ефекту аліасингу, використовується спеціальний фільтр згладжування (фільтр низьких частот) встановлюється перед АЦП. Він пропускає частоти, нижчі за половину частоти дискретизації АЦП, і відсікає вищі частоти.
Щоб обчислити спектр сигналу за його дискретними відліками, дискретний Дискретне перетворення Фур'є (DFT) використовується. Зауважимо ще раз, що спектр дискретного сигналу “за визначенням” обмежений частотою Fmax, меншою за половину частоти дискретизації Fd. Тому спектр дискретного сигналу можна подати як суму скінченний кількість гармонік, на відміну від нескінченної суми для ряду Фур'є неперервного сигналу, спектр якого може бути необмеженим. Згідно з теоремою Котельникова, максимальна частота гармоніки повинна бути такою, щоб на неї припадало не менше двох відліків, тому кількість гармонік дорівнює половині кількості відліків дискретного сигналу. Тобто, якщо у вибірці N відліків, то кількість гармонік у спектрі буде N/2.
Розглянемо тепер дискретне перетворення Фур'є (ДПФ).
Порівняння з рядом Фур'є
Як бачимо, вони збігаються, за винятком того, що час у ШПФ дискретний, а кількість гармонік обмежена N/2, що вдвічі менше за кількість відліків.
Формули ДПФ записуються в безрозмірних цілочисельних змінних k, s, де k - кількість відліків сигналу, s - кількість спектральних компонент.
Величина s показує кількість повних гармонічних коливань за період T (тривалість вимірювання сигналу). Дискретне перетворення Фур'є використовується для знаходження амплітуд і фаз гармонік чисельно, тобто "на комп'ютері".
Як вже було сказано вище, при розкладанні неперіодичної функції (нашого сигналу) в ряд Фур'є, отриманий ряд Фур'є фактично відповідає періодичній функції з періодом T (рис.12).
Рис.12. Періодична функція f(x) з періодом T0, з періодом T>T0
Як видно з рис. 12, функція f(x) є періодичною з періодом T0. Однак через те, що довжина вимірювального відрізка T не дорівнює періоду функції T0, функція, отримана у вигляді ряду Фур’є, має розрив у точці T. Як наслідок, спектр цієї функції міститиме велику кількість високочастотних гармонік. Це явище відоме як спектральний витік, а на практиці це скорочується до віконування сигнал до перетворення. Якщо тривалість вимірювального відриву T збігалася з періодом функції T0, то спектр, отриманий після перетворення Фур’є, містив би лише першу гармоніку (синусоїду з періодом, рівним тривалості відриву), оскільки функція f(x) є синусоїдою.
Іншими словами, програма ДПФ “не знає”, що наш сигнал є “фрагментом синусоїди”, а намагається подати у вигляді ряду періодичну функцію, яка має розрив через неузгодженість окремих фрагментів синусоїди.
В результаті в спектрі з'являються гармоніки, які в сукупності повинні відображати форму функції, включаючи цей розрив.
Таким чином, щоб отримати “правильний” спектр сигналу, який є сумою кількох синусоїд з різними періодами, необхідно, щоб ціла кількість періодів кожна синусоїда повинна бути присутня на періоді вимірювання сигналу. На практиці ця умова може бути виконана при достатньо великій тривалості вимірювання сигналу.
Рис.13 Приклад функції та спектру сигналу кінематичної похибки редуктора
За меншої тривалості картина виглядатиме “гірше”:
Рис.14 Приклад функції та спектру вібрації ротора
На практиці буває важко зрозуміти, де “реальні складові”, а де “артефакти”, спричинені неузгодженістю періодів складових і тривалості дискретизації сигналу або “стрибками та розривами” у формі сигналу. Звісно, слова “реальні складові” і “артефакти” взяті в лапки не випадково. Наявність багатьох гармонік на графіку спектра не означає, що наш сигнал справді складається з них. Це все одно що думати, ніби число 7 “складається” з чисел 3 і 4. Число 7 можна подати як суму 3 і 4 – і це правильно.
Так само і наш сигнал… точніше, навіть не “наш сигнал”, а періодична функція, утворена повторенням нашого сигналу (вибірки), може бути подана як сума гармонік (синусоїд) з певними амплітудами і фазами. Але в багатьох важливих для практики випадках (див. рисунки вище) гармоніки, отримані у спектрі, справді можна пов'язати з реальними процесами, що мають циклічний характер і суттєво впливають на форму сигналу.
Деякі результати
1. Реальний вимірюваний сигнал тривалістю T секунд, оцифрований АЦП, тобто представлений набором дискретних відліків (N штук), має дискретний неперіодичний спектр, представлений набором гармонік (N/2 штук).
2. Сигнал представлений набором дійсних значень. Його спектр DFT є набором комплексних коефіцієнтів зі спряженою симетрією; з них отримують амплітудний спектр — набір дійсних невід’ємних амплітуд (і фаз) на додатних частотах, і саме цей однобічний амплітудний спектр на практиці відображається на графіку. Двобічна комплексна форма з від’ємними частотами та однобічна форма амплітуди/фази є еквівалентними представленнями того самого спектра — для аналізу сигналів зазвичай зручніше працювати з однобічним амплітудним спектром.
3. Сигнал, виміряний протягом часу T, визначений лише на цьому інтервалі T. Що відбувалося до початку вимірювання і що буде після нього, науці невідомо, та в нашому випадку це нецікаво. ШПФ обмеженого в часі сигналу дає його “реальний” спектр у тому сенсі, що за певних умов дозволяє обчислити амплітуди й частоти його складових.
