Разбираемся в ДПФ

Дискретное преобразование Фурье - краеугольный камень цифровой обработки сигналов. Формула ДПФ выглядит следующим образом:

N - количество отсчетов ДПФ.

m = 0..N-1 - индекс ДПФ в частотной области.

 n = 0..N-1 - временной индекс входных отсчетов.

x(n) - входные отсчеты во временной области.

X(m) - выходные отсчеты в частотной области.

Начнем, пожалуй, с самого сложного слагаемого в этой сумме:

По формуле Эйлера комплексную экспоненту можно представить как сумму тригонометрических функций. Рассмотрим как меняются индексы m и n при ДПФ.

Рассмотрим пример для N = 4. Формула ДПФ тогда выглядит следующим образом:

Почленно:

X(0) = x(0)*(cos(2*Pi*0*0/4) - jsin(2*Pi*0*0/4))

        + x(1)*(cos(2*Pi*1*0/4) - jsin(2*Pi*1*0/4))

        + x(2)*(cos(2*Pi*2*0/4) - jsin(2*Pi*2*0/4))

        + x(3)*(cos(2*Pi*3*0/4) - jsin(2*Pi*3*0/4)) 

        =

           x(0)*(cos(0*Pi) - jsin(0*Pi))

        + x(1)*(cos(0*Pi) - jsin(0*Pi))

        + x(2)*(cos(0*Pi) - jsin(0*Pi))

        + x(3)*(cos(0*Pi) - jsin(0*Pi))

X(1) = x(0)*(cos(2*Pi*0*1/4) - jsin(2*Pi*0*1/4))

        + x(1)*(cos(2*Pi*1*1/4) - jsin(2*Pi*1*1/4))

        + x(2)*(cos(2*Pi*2*1/4) - jsin(2*Pi*2*1/4))

        + x(3)*(cos(2*Pi*3*1/4) - jsin(2*Pi*3*1/4))

        =

           x(0)*(cos(0*Pi) - jsin(0*Pi))

        + x(1)*(cos(Pi/2) - jsin(Pi/2))

        + x(2)*(cos(Pi) - jsin(Pi))

        + x(3)*(cos(3*Pi/2) - jsin(3*Pi/2))

X(2) = x(0)*(cos(2*Pi*0*2/4) - jsin(2*Pi*0*2/4))

        + x(1)*(cos(2*Pi*1*2/4) - jsin(2*Pi*1*2/4))

        + x(2)*(cos(2*Pi*2*2/4) - jsin(2*Pi*2*2/4))

        + x(3)*(cos(2*Pi*3*2/4) - jsin(2*Pi*3*2/4))

        =

           x(0)*(cos(0*Pi) - jsin(0*Pi))

        + x(1)*(cos(Pi) - jsin(Pi))

        + x(2)*(cos(2*Pi) - jsin(2*Pi))

        + x(3)*(cos(3*Pi) - jsin(3*Pi))

X(3) = x(0)*(cos(2*Pi*0*3/4) - jsin(2*Pi*0*3/4))

        + x(1)*(cos(2*Pi*1*3/4) - jsin(2*Pi*1*3/4))

        + x(2)*(cos(2*Pi*2*3/4) - jsin(2*Pi*2*3/4))

        + x(3)*(cos(2*Pi*3*3/4) - jsin(2*Pi*3*3/4))

        =

           x(0)*(cos(0*Pi) - jsin(0*Pi))

        + x(1)*(cos(3*Pi/2) - jsin(3*Pi/2))

        + x(2)*(cos(3*Pi) - jsin(3*Pi))

        + x(3)*(cos(9*Pi/2) - jsin(9*Pi/2))

Можно заметить, что всегда с входным отсчетом x(0) аргумент комплексной экспоненты равен нулю. При х(1) аргумент равен m*Pi/2. При x(2) - m*Pi, x(3) - m*3*Pi/2. Другими словами: для каждого m-ного выходного отсчета аргумент у тригонометрических функций совершает m - оборотов, причем прирост аргумента величина постоянная и равная 2*Pi*m/N. Какой в этом смысл? Это фактически перенос на нулевую частоту, можно посмотреть здесь как это выводится. Номер выходного отсчета связан с частотой ДПФ выражением:

fs - частота дискретизации сигнала. Так для случая ДПФ4 при частоте дискретизации 4 кГц, анализируемые частоты в спектре будут для X(0) - 0 кГц, X(1) - 1, X(2) - 2, X(3) - 3. Соответственно если смотреть на аргумент у выходных отсчетов, то при X(1) аргумент делает 1 оборот за 4 отсчета, что соответствует 1 кГц. Если не очень понятно, то можно почитать тут подробнее.
Поскольку функция синуса и косинуса периодическая, то для любого аргумента справедливо cos(2*Pi+m) = cos(m). Так же будем пользоваться четностью косинуса и нечетностью синуса. Упростить выражения для выходных отсчетов:

X(0) = x(0) + x(1) + x(2) + x(3) = (x(0) + x(2)) + (x(1) + x(3))

X(1) = x(0)

        + x(1)*(cos(Pi/2) - jsin(Pi/2))

        - x(2)

        + x(3)*(cos(3*Pi/2) - jsin(3*Pi/2))

        =

           x(0) - x(2) + jx(1) -jx(3) = (x(0)-x(2)) - j(x(1)-x(3))

X(2) = x(0)

        - x(1)

        + x(2)

        - x(3)

        =

         (x(0)+x(2))-(x(1)+x(3))

X(3) =  x(0)

        + x(1)*(cos(3*Pi/2) - jsin(3*Pi/2))

        -  x(2)

        + x(3)*(cos(Pi/2) - jsin(Pi/2))

        =

           (x(0) - x(2)) + j(x(1)-x(3)))

Теперь нужно проделать те же действия для N = 8. Получится следующий результат:

Отметим, что входная последовательность с АЦП всегда действительная, поэтому ДПФ симметричен относительно середины: m отсчет ДПФ имеет тот же модуль, что и N-m. Фазовые углы m отсчета равен N-m взятым со знаком минус. Таким образом получается, что пара выходных отсчетов ДПФ m и N-m комплексно сопряжены, но только в случае действительной входной последовательности. Это означает, что в N точечном ДПФ только N/2+1 отсчетов независимы и для вычисления ДПФ достаточно вычислить только их.
Вот простенькая реализация ДПФ на си:

#include <stdio.h> #include <math.h> int const N = 8; void DFT(double * pIN, double * pOUT, int N) { double delta_f = 2*M_PI/N; double RE_out[N]; double IM_out[N]; for(int i = 0; i < N; i++) { RE_out[i] = .0; IM_out[i] = .0; } for(int i = 0; i < N; i++) { for(int j = 0; j < N; j++) { RE_out[i] += *(pIN+j)*cos(delta_f*i*j); IM_out[i] += *(pIN+j)*sin(delta_f*i*j)* -1.0; } *(pOUT+i) = sqrt(RE_out[i]*RE_out[i] + IM_out[i]*IM_out[i]); printf("%d\t\tRE: %5.1f\t IM: %5.1f\tAMPL: %5.1f\n",i, RE_out[i], IM_out[i], *(pOUT+i)); } } int main(){ double * pIN = new double[N]; double * pOUT = new double[N]; for(int i = 0; i < N; i++) { *(pIN+i) = cos(2*M_PI*1*i/N); *(pOUT+i) = .0; } DFT(pIN, pOUT, N); return 0; }

Дальше будем разбираться как посчитать БПФ