Алгоритмы цифровой обработки сигналов (ЦОС): БПФ (FFT), цифровые фильтры (реализация на C)

Введение

Цифровая обработка сигналов (ЦОС, DSP — Digital Signal Processing) — основа современных систем: от аудио/видео до радиолокации и биомедицины. Ключевые задачи:

• анализ спектра сигнала (БПФ);
• фильтрация помех и выделение полезных компонент;
• сжатие и кодирование данных;
• обнаружение и распознавание сигналов.

В статье разберём:

• основы дискретизации и теоремы Котельникова;
• алгоритм быстрого преобразования Фурье (FFT);
• типы цифровых фильтров (FIR, IIR);
• практическую реализацию на языке C для микроконтроллеров и ПК.

---

1. Основы дискретизации и представления сигналов

1.1. Теорема Котельникова (Найквиста‑Шеннона)

Чтобы аналоговый сигнал можно было точно восстановить по дискретным отсчётам:

• частота дискретизации fs​ должна быть не менее чем вдвое выше максимальной частоты сигнала fmax​:fs​≥2⋅fmax​
• иначе возникает алиасинг (наложение спектров).

Пример: для аудио до 20 кГц нужна fs​≥40 кГц (стандарт — 44.1 кГц или 48 кГц).

1.2. Представление сигнала в ЦОС

• Дискретный сигнал: последовательность отсчётов x[n], где n=0,1,2,…,N−1.
• Длина блока: N отсчётов (для БПФ часто N=2k).
• Разрешение по частоте после БПФ:Δf=Nfs​​

1.3. Квантование и шум

• АЦП преобразует аналоговое напряжение в целое число (например, 12‑битное).
• Погрешность квантования → шум (обычно моделируется как белый шум).
• Для снижения шума применяют передискретизацию (oversampling).

---

2. Быстрое преобразование Фурье (FFT)

2.1. Суть и назначение

FFT — эффективный алгоритм вычисления дискретного преобразования Фурье (ДПФ):

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

ДПФ (прямое):

X[k]=n=0∑N−1​x[n]⋅e−jN2π​kn,k=0,1,…,N−1

FFT сокращает сложность с O(N2) до O(Nlog2​N).

2.2. Алгоритм Кули‑Тьюки (Radix‑2)

Основные шаги:

1. Перестановка отсчётов (bit‑reversal): индексы переставляются в «зеркальном» порядке.
2. Бабочка (butterfly operation): базовая операция для пары точек:{A′=A+B⋅WB′=A−B⋅W​ где W — комплексный поворот (twiddle factor): W=e−jN2π​k.
3. Итерации по уровням: на каждом уровне обрабатывается блок размером 2m.

2.3. Реализация FFT на C (упрощённая версия)

#include <math.h>
#include <complex.h>

#define PI 3.14159265358979323846


// Комплексное число (реальная и мнимая часть)
typedef struct {
    double re;
    double im;
} complex_t;

// Бабочка FFT
void butterfly(complex_t *a, complex_t *b, complex_t w) {
    complex_t t;
    t.re = w.re * b->re - w.im * b->im;
    t.im = w.re * b->im + w.im * b->re;

    b->re = a->re - t.re;
    b->im = a->im - t.im;
    a->re += t.re;
    a->im += t.im;
}

// FFT (Radix-2, N должно быть степенью 2)
void fft(complex_t *x, int N) {
    // 1. Bit-reversal перестановка
    for (int i = 0; i < N; i++) {
        int j = 0;
        for (int k = 0; k < log2(N); k++) {
            j = (j << 1) | ((i >> k) & 1);
        }
        if (j > i) {
            complex_t tmp = x[i];
            x[i] = x[j];
            x[j] = tmp;
        }
    }

    // 2. Итерации по уровням
    for (int s = 1; s <= log2(N); s++) {
        int m = 1 << s;            // Размер блока
        int half_m = m >> 1;
        complex_t w;              // Поворотный множитель
        w.re = cos(-2.0 * PI / m);
        w.im = sin(-2.0 * PI / m);

        for (int k = 0; k < N; k += m) {
            complex_t w_power = {1.0, 0.0};  // w^0
            for (int j = 0; j < half_m; j++) {
                butterfly(&x[k + j], &x[k + j + half_m], w_power);
                // Обновляем w_power = w^(j+1)
                complex_t tmp;
                tmp.re = w_power.re * w.re - w_power.im * w.im;
                tmp.im = w_power.re * w.im + w_power.im * w.re;
                w_power = tmp;
            }
        }
    }
}

2.4. Практические аспекты

• Окно (windowing): перед FFT применяют окна (Хэмминга, Ханна) для снижения спектральных утечек.
• Нормализация: амплитуды после FFT нужно делить на N/2 (для вещественного сигнала).
• Обратное FFT (IFFT): аналогично, но с положительным знаком в экспоненте и делением на N.

---

3. Цифровые фильтры

3.1. Типы фильтров

1. FIR (Finite Impulse Response) — с конечной импульсной характеристикой:
   • линейная фаза;
   • устойчивость (всегда);
   • высокая вычислительная сложность.
2. IIR (Infinite Impulse Response) — с бесконечной импульсной характеристикой:
   • нелинейная фаза;
   • возможная неустойчивость;
   • низкая сложность (меньше коэффициентов).

3.2. FIR‑фильтры

Уравнение:

y[n]=k=0∑M​bk​⋅x[n−k]

где:

• x[n] — входной сигнал;
• y[n] — выходной сигнал;
• bk​ — коэффициенты фильтра;
• M — порядок фильтра.

Проектирование:

• метод окон (Хэмминга, Блэкмана);
• метод частотной выборки;
• оптимизация (алгоритм Паркса‑Макклеллана).

Реализация на C:

#define FIR_ORDER 64

double fir_filter(double x, double *coeffs, double *delay_line, int *idx) {
    delay_line[*idx] = x;  // Сохраняем новый отсчет
    double y = 0.0;
    
    // Свёртка
    for (int i = 0; i <= FIR_ORDER; i++) {
        y += coeffs[i] * delay_line[(*idx - i + FIR_ORDER + 1) % (FIR_ORDER + 1)];
    }
    
    *idx = (*idx + 1) % (FIR_ORDER + 1);  // Сдвигаем индекс
    return y;
}

3.3. IIR‑фильтры

**