Розрахунок і реалізація цифрових фільтрів. Просто і швидко

Корисна інформація для тих, хто програмно реалізовує цифрову фільтрацію сигналу.
Програмна реалізація цифрових фільтрів, алгоритми і приклади підпрограм цифрової фільтрації.

Цифрова обробка сигналів, зокрема, цифрова фільтрація сигналів, є частою і популярною задачею для інженерів. Якщо потрібно розрахувати необхідний цифровий фільтр, то в наш час не варто робити це самостійно, з допомогою калькулятора або самостійно написаної програми. Досить скористатись для цього давно розробленим програмним інструментом.
Ви не зробите революцію в проектуванні цифрових фільтрів, навіть не надійтесь. Люди з прізвищами Баттерворт, Знов’як, Чебишев та інші вже детально обкатали цю тему, а вам краще скористатись готовими результатами.

Спочатку коротко розглянемо розрахунок цифрових фільтрів з допомогою популярної програми MATLAB. Потім розглянемо практичне використання розрахованих коефіцієнтів потрібного нам фільтра, тобто програмні алгоритми для цифрової фільтрації сигнала.
Хтось скаже, що MATLAB є платною програмою. Цілком можливо, що навіть платною, хоча завдяли Інтернету є можливість цього не відчувати.

Що робить підпрограма (любителі C# називають “метод”), яка є моделлю цифрового фільтра ? За кожним звертанням до підпрограми їй передається одне поточне вхідне значення параметра, який фільтруємо, а підпрограма повертає вирахуване поточне вихідне значення. Так само, як апаратний фільтр: щось нескінченно подаємо на вхід, і відфільтроване знімаємо з вихода.

MATLAB пропонує підтримку для цілочисельної математики і математики з плаваючою точкою, а також мовні конструкції для обробки і аналізу великих наборів даних. Велика кількість оптимізацій для типів даних, операцій, функцій і апаратного забезпечення привело до значного підвищення швидкості обчислень. Дуже корисним є те, що MATLAB чудово вираховує коефіцієнти цифрового фільтра, які потім використовуються в імітуючій підпрограмі.

Серед програмних можливостей MATLAB слід зазначити зручний англомовний інтерфейс, що включає нові програмні і налагоджувальні інструменти, автоматичний аналіз якості коду, а також можливість збереження файлу програми безпосередньо в HTML- або Word- форматах. Нові інтерактивні інструменти побудови графіків, включаючи можливість генерації коду для повторного багатократного створення графіків.

Цілочисельна математика і математика одинарної точності. Обробка нових типів даних без перетворення їх в числа подвійної точності значно підвищує продуктивність і зменшує об’єм використовуваної пам’яті. Це дає можливість працювати з великими масивами даних. За допомогою FFT-алгоритмів (Fast Fourier Transform, швидке перетворення Фур’є) швидкість  Фур’є-перетворень одинарної точності підвищена в середньому на 20%. окрім того, MATLAB використовує бібліотеку цілочисельних алгоритмів Intel MMX, що підвищує швидкість обчислень з цілочисельними даними до 8 разів.

Компілятор MATLAB Compiler підтримує усю мову MATLAB, включаючи більшість додатків (MATLAB Toolboxes).
Поліпшений компілятор MATLAB Compiler дає можливість інженерам поширювати незалежні застосування, розроблені в MATLAB, або включати їх в такі засоби розробки, як Excel, C  і COM. В результаті, інженери і вчені тепер зможуть творити значно ширший спектр MATLAB-додатків і їх поширювати.

Що міститься в сімействі продуктів MATLAB ?
MATLAB
Simulink
Aerospace Blockset
Bioinformatics Toolbox
CDMA Reference Blockset
— та багато іншого, це залежить від версії MATLAB.

Отже, розраховуємо фільтр.

Запустили MATLAB. В вікні Command Window потрібно ввести команду “fdatool”, після чого запуститься окрема програма для розрахунку фільтрів.
Вид фільтра вибираємо:
LOWPASS filter (LPF) – фільтр низької частоти (ФНЧ)
HIGHPASS filter (HPF) – фільтр високої частоти (ФВЧ)
BANDPASS filter (BPF) – смуговий фільтр
BANDSTOP filter (BSF) – загороджувальний фільтр

Також можна вибрати такі математичні моделі для цифрового перетворення сигналу:
DIFFERENTIATOR – диференціатор
MULTIBAND – багатосмуговий фільтр
HILBERT TRANSFORMER – перетворювач Гільберта. амплітуда сигнала не змінюється, але зміщується фаза кожної синусоїди на ± 90 °.
ARBITRARY MAGNITUDE – довільна АЧХ (амплітудно-частотна характеристика)
ARBITRARY GROUP DELAY – довільний фільтр групової затримки, такі фільтри можна використовувати для корекції фазових спотворень, введених іншими фільтрами.
PEEKING – вузькополосий фільтр. Застосовується, щоб зберегти певну частотну полосу
NOTCHING – режекторний фільтр. Застосовується, щоб усунути певну частотну полосу

Можна вибирати фільтр з кінцевою імпульсною рарактеристикою FIR (Finite impulse response), або з нескінченною імпульсною рарактеристикою IIR (Infinite impulse response). Фільтри IIR традиційно дають кращі результати, саме їх будемо розглядати. Встановлюємо параметри фільтра, вибираємо метод розрахунку (Баттерворта, Чебишева, та інші), порядок фільтра можна вибирати самому або вказати пункт “Мінімізувати порядок”  (Minimum order).

Розглянемо приклад. Вважаємо, що ми розрахували з допомогою MATLAB цифровий фільтр типу IIR. Як одержати необхідні коефіцієнти фільтра ?
Можна через Menu – File – Export – Coeficient File одержати файл “untitled.fcf” в якому ми побачимо, наприклад, такі числа:

———————————————————————————-
Generated by MATLAB(R) and the Signal Processing Toolbox

Generated on:

Discrete-Time IIR Filter (real)
——————————-
Filter Structure  : Direct-Form II, Second-Order Sections
Filter Order      : 20
Stable            : Yes
Linear Phase      : No

SOS matrix:
1  0  -1  1  -1.945902111935361  0.974222744088694
1  0  -1  1  -1.987674236609505  0.991043181597264
1  0  -1  1  -1.902652268762338  0.928321966443545
1  0  -1  1  -1.969177121612932  0.972775889099217
1  0  -1  1  -1.873487014461808  0.895037788778907
1  0  -1  1  -1.948824410675406  0.952998193881200
1  0  -1  1  -1.925123116888261  0.930458301381064
1  0  -1  1  -1.862815343397140  0.879358495490341
1  0  -1  1  -1.897690896486521  0.905289876215382
1  0  -1  1  -1.872688833141580  0.884190212397096
Scale Factors:
0.055709918470064
0.055709918470064
0.054817769783290
0.054817769783290
0.054091680572213
0.054091680572213
0.053581879575375
0.053581879575375
0.053319577182277
0.053319577182277
1.000000000000000

В кожному рядку масива  SOS matrix ми бачимо 6 чисел. В нашому прикладі вийшло 10 рядків, а в таблиці Scale Factors  – 11 рядків. Формуємо з таблиці SOS matrix двомірний масив   A[MWSPT_NSEC][6].
Параметр MWSPT_NSEC – це кількість горизонтальних рядків в матриці SOS matrix, в нашому випадку = 10.
Ось так для компілятора “C” чи “C++”  виглядатиме цей двохмірний масив:

#define MWSPT_NSEC 10        // тому що 10 рядків

// якщо користуємось C#, нам недоступна #define для такого використання,
// тому замість MWSPT_NSEC використовуємо числове значення, в даному випадку 10

// зформований масив із матриці SOS matrix:
const float A[MWSPT_NSEC][6] = {
1, 0, -1, 1,  -1.945902111935361,  0.974222744088694,
1, 0, -1, 1,  -1.987674236609505,  0.991043181597264,
1, 0, -1, 1,  -1.902652268762338,  0.928321966443545,
1, 0, -1, 1,  -1.969177121612932,  0.972775889099217,
1, 0, -1, 1,  -1.873487014461808,  0.895037788778907,
1, 0, -1, 1,  -1.948824410675406,  0.952998193881200,
1, 0, -1, 1,  -1.925123116888261,  0.930458301381064,
1, 0, -1, 1,  -1.862815343397140,  0.879358495490341,
1, 0, -1, 1,  -1.897690896486521,  0.905289876215382,
1, 0, -1, 1,  -1.872688833141580,  0.884190212397096 };

З масиву “Scale Factors” формуємо одномірний масив B[MWSPT_NSEC].

const float B[MWSPT_NSEC] ={
0.055709918470064,
0.055709918470064,
0.054817769783290,
0.054817769783290,
0.054091680572213,
0.054091680572213,
0.053581879575375,
0.053581879575375,
0.053319577182277,
0.053319577182277};

А останнє число в  масиві “Scale Factors”  – це перемножуючий коефіцієнт C.
const float C =  1.000000000000000 ;
Якщо MATLAB не видає це останнє число, а кількість рідків в SOS matrix дорівнює кількості рядків в Scale Factors, значить, const float C =  1.0

Для програмної реалізації фільтра можна скористатись спрощеним алгоритмом Знов’яка (Znoviak algorithm).

Алгоритм не є абсолютно універсальним, але він швидкий, і чудово підходить для фільтрів низької частоти, високої частоти, смугового та загороджувального фільтра.
При цьому підпрограма, написана на “C“, виглядає так:

float IIR_ZNO_filter (float X)            // Znoviak algorithm
{
#define MWSPT_NSEC 10

const float A[MWSPT_NSEC][6] = {
1, 0, -1, 1,  -1.945902111935361,  0.974222744088694,
1, 0, -1, 1,  -1.987674236609505,  0.991043181597264,
1, 0, -1, 1,  -1.902652268762338,  0.928321966443545,
1, 0, -1, 1,  -1.969177121612932,  0.972775889099217,
1, 0, -1, 1,  -1.873487014461808,  0.895037788778907,
1, 0, -1, 1,  -1.948824410675406,  0.952998193881200,
1, 0, -1, 1,  -1.925123116888261,  0.930458301381064,
1, 0, -1, 1,  -1.862815343397140,  0.879358495490341,
1, 0, -1, 1,  -1.897690896486521,  0.905289876215382,
1, 0, -1, 1,  -1.872688833141580,  0.884190212397096 };

const float B[MWSPT_NSEC] ={
0.055709918470064,
0.055709918470064,
0.054817769783290,
0.054817769783290,
0.054091680572213,
0.054091680572213,
0.053581879575375,
0.053581879575375,
0.053319577182277,
0.053319577182277 };

const float C =  1.000000000000000 ;

static float W[MWSPT_NSEC][3];
//
// це внутрішній масив підпрограми, обов’язково “static”
//
// пояснення для тих, хто не любить програмувати на “C“:
// признак “static” означає, що значення масиву W повинні зберігатись
// при виході з підпрограми, не бути тимчасовими,
// або навіть регістровими. Саме тому при використанні C# масив
// розміщується за межами “метода” (підпрограми)
//
for (int i = 0; i < MWSPT_NSEC; i++)
{
W[i][2] = W[i][1];
W[i][1] = W[i][0];
W[i][0] = X * B[i] – W[i][1] * A[i][4] – W[i][2] * A[i][5];
X = W[i][0] * A[i][0] + W[i][1] * A[i][1] + W[i][2] * A[i][2];
}
X=X*С;
return X;
}

При кожному звертанні до підпрограми “IIR_ZNO_filter” передається ПОТОЧНЕ значення вхідного дискретизованого сигнала X, а на виході підпрограма видає відфільтроване поточне значення X.
Зрозуміло, що, так само як і в апаратному фільтрі, спочатку йде перехідний процес, тому перші відліки відфільтрованого сигнала можуть бути некоректні.

А тепер, подивившись на коефіцієнти нашого фільтра, можемо застосувати хитрість і зменшити кількість обчислень..Подивіться на матрицю А, в нашому конкретному прикладі цифрового фільтра кожен рядок починається числами 1, 0, -1.
Навіщо нам перемножувати що-небудь на  1, або 0, або -1 ? Обійдемось, скоротимо алгоритм конкретно для нашого фільтра. Змінимо лише один рядок.

 // X = W[i][0] * A[i][0] + W[i][1] * A[i][1] + W[i][2] * A[i][2];  // так було
X = W[i][0] * A[i][0] – W[i][2] * A[i][2];                                    // так стало

Це лише один із способів реалізації цифрового фільтра. Можна піти іншим шляхом.
В програмі FDATool через Menu – Targets – Generate C Header можна одержати готовий файл з коефіцієнтами під назвою “fdacoefs.h”.
Ось так виглядає приклад цього файла:

/*
* Discrete-Time IIR Filter (real)
* ——————————-
* Filter Structure  : Direct-Form II, Second-Order Sections
* Filter Order      : 20
* Stable              : Yes
* Linear Phase   : No
*/

/* General type conversion for MATLAB generated C-code  */
#include “tmwtypes.h”

#define MWSPT_NSEC 21
const int NL[MWSPT_NSEC] = { 1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1 };
const float NUM[MWSPT_NSEC][3] = {
{    0.05570991847006,                 0,                 0   },
{                   1,                 0,                -1  },
{    0.05570991847006,                 0,                 0   },
{                   1,                 0,                -1   },
{    0.05481776978329,                 0,                 0   },
{                   1,                 0,                -1   },
{    0.05481776978329,                 0,                 0   },
{                   1,                 0,                -1   },
{    0.05409168057221,                 0,                 0   },
{                   1,                 0,                -1   },
{    0.05409168057221,                 0,                 0   },
{                   1,                 0,                -1  },
{    0.05358187957537,                 0,                 0   },
{                   1,                 0,                -1   },
{    0.05358187957537,                 0,                 0   },
{                   1,                 0,                -1  },
{    0.05331957718228,                 0,                 0   },
{                   1,                 0,                -1   },
{    0.05331957718228,                 0,                 0   },
{                   1,                 0,                -1   },
{                   1,                 0,                 0   }
};
const int DL[MWSPT_NSEC] = { 1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1 };
const float DEN[MWSPT_NSEC][3] = {
{                   1,                 0,                 0   },
{                   1,   -1.945902111935,   0.9742227440887   },
{                   1,                 0,                 0   },
{                   1,    -1.98767423661,   0.9910431815973   },
{                   1,                 0,                 0   },
{                   1,   -1.902652268762,   0.9283219664435   },
{                   1,                 0,                 0   },
{                   1,   -1.969177121613,   0.9727758890992   },
{                   1,                 0,                 0   },
{                   1,   -1.873487014462,   0.8950377887789   },
{                   1,                 0,                 0   },
{                   1,   -1.948824410675,   0.9529981938812   },
{                   1,                 0,                 0   },
{                   1,   -1.925123116888,   0.9304583013811  },
{                   1,                 0,                 0   },
{                   1,   -1.862815343397,   0.8793584954903   },
{                   1,                 0,                 0   },
{                   1,   -1.897690896487,   0.9052898762154   },
{                   1,                 0,                 0   },
{                   1,   -1.872688833142,   0.8841902123971   },
{                   1,                 0,                 0   }
};

Підпрограма, яка використовує цей файл, виглядає так:

float IIR_filt (float x)
{
#include “fdacoefs.h”

int i,j;
float w;
static float WHistory[MWSPT_NSEC][2]; // це внутрішній масив, обов’язково static
// (при використанні C# цей масив повинен бути за межами “метода”, так прийнято
// називати підпрогаму, якщо користуємось C#)

for ( i = 0; i < MWSPT_NSEC; i++){
w = x;
for ( j=1; j<DL[i]; j++){w-= DEN[i][j] * WHistory[i][j-1];};
w*=1/DEN[i][0];

x = w * NUM[i][0];
for ( j=1; j<NL[i]; j++){x+=NUM[i][j] * WHistory[i][j-1];};
WHistory[i][1] = WHistory[i][0];
WHistory[i][0] = w;
};
return x;
}

Так само при кожному звертанні до підпрограми передається ПОТОЧНЕ значення вхідного дискретизованого сигнала x, а на виході підпрограма видає відфільтроване поточне значення x. Цей алгоритм виконує трохи більше обчислень для виконання тої самої задачі, і в деяких випадках він є більш універсальним.
Залишився невідомим файл tmwtypes.h але цей файл буде завантажений, якщо натиснути на його назву.

Нічого складного. Щось незрозуміло – запитайте в коментарях чи на форумі, я поясню, заодно доповню статтю.

Весь список статей нашого сайту, корисних для програмістів, знаходиться ось тут