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

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

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; }

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

Источник синусоидального сигнала

Продолжаем искать применение для CORDIC. Уже была демонстрация различных видов модуляции, теперь сделаем программируемый источник синусоидального сигнала. Программируемым он будет потому, что можно задать частоту выходного сигнала. Генерировать будем и синус, и косинус.

Нам нужно воспроизвести прирост аргумента, равный 2*Pi*f/F, где f - частота выходного сигнала; F - частота дискретизации. Частоту дискретизации возьмем равной 8 кГц, что соответствует одному речевому каналу. Константа 2*Pi для CORDIC была посчитана, операция деления реализована, частота f задается извне. Еще нам понадобятся временные метки, в нашем случае идущие с частотой 8 кГц. Я воспользуюсь бодовым генератором.

Проекты со временем разрастаются, поэтому постепенно я начну переезд на гитхаб. Ну что же, все основные ссылки даны, теперь картинка для привлечения внимания:

Меняем значение на входе out_freq, меняется частота на выходах sin и cos

Какие неявные ограничения есть в такой реализации? Во-первых частота дискретизации равна 8 кГц, это означает, что мы не можем генерировать частоту выше 4 кГц. Частоту повышать будем с помощью интерполирующих фильтров, но об этом потом.

Список интересующих файлов с гитхаба:

cordic_gen.v  - топовый модуль, реализующий логику приращения аргумента для модуля CORDIC.

divider.v - модуль деления с остатком, где делимое 2*Pi*f, делитель F. Смысл остатка для топового модуля следующий: как только счетчик досчитывает до значения остатка от деления, приращение, равное частному, инкрементируется на единицу и счетчик сбрасывается. 

main.v - модуль, реализующий CORDIC. 

Деление на Verilog

Я уже рассматривал операции сложения, вычитания и умножения. Так вышло, что операция деления не является базовой в языке Verilog, поэтому нельзя делить одну регистровую переменную на другую. Сегодня я попытаюсь понять, почему так произошло, и как делить в целочисленных.
Первым делом я начал гуглить этот вопрос и наткнулся на довольно интересную страницу, где описывается умножение и деление на Verilog. Сразу понимания сути это описание не принесло, поэтому стал искать книгу из списка литературы "Computer Organization & Design". И вот в этой книге уже довольно подробно расписано что и почему. Очень советую полистать книгу на досуге.
Сначала определим с чем мы будем иметь дело по ходу работы. У нас будут делимое (divident), делитель (divider), частное (quotient) и остаток (reminder). Причем делимое = делитель*частное + остаток. Для простоты будем рассматривать случай беззнаковых 16 битных целых чисел. Алгоритм выглядит следующим образом:
0. Инициализируем частное = 0. Поскольку каждую итерацию алгоритма делитель нужно будет сдвигать на один разряд вправо, то сделаем 32 битную "копию" делителя, добавив в младшие разряды 16 нулей. Остаток приравниваем к делимому. Дальше весь алгоритм завязан на частном, делителе и остатке. Алгоритм займет N+1 итерацию, где N - количество разрядов делимого и в нашем случае равно шестнадцати.
1. Остаток = остаток - делитель. Если остаток меньше нуля, то вернуть прежний остаток (остаток = остаток + делитель), сдвинуть частное влево на один разряд и в младший разряд подставить 0.
2. Сдвинуть делитель вправо на один разряд. Повторить п. 1.
3..17. Посторить п.2.

Рассмотрим пример из книжки. Поделим 8 битное число 7 на 4 битное число 2.
divident = 7 = 0000 0111;
divider = 2 = 0010;

0. Инициализация:
quotient = 0000;
divider = 0010 0000;
reminder = 0000 0111;

1.
reminder = reminder - divider = 1110 0111 < 0;
reminder = 0000 0111;
quotient = quotient << 1 + 0 = 0000;

2.
divider = divider >> 1 = 0001 0000;
reminder = reminder - divider = 11110111 < 0;
reminder = 0000 0111;
quotient = quotient << 1 + 0 = 0000;

3.
divider = divider >> 1 = 0000 1000;
reminder = reminder - divider = 11111111 < 0;
reminder = 0000 0111;
quotient = quotient << 1 + 0 = 0000;

4.
divider = divider >> 1 = 0000 0100;
reminder = reminder - divider = 0000 0011 > 0;
quotient = quotient << 1 +1 = 0001;

5.
divider = divider >> 1 = 0000 0010;
reminder = reminder - divider = 0000 0001 > 0;
quotient = quotient << 1 +1 = 0011;
end

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

Ну и теперь реализация на Verilog. Количество итераций здесь равно разрядности делимого и получено это благодаря тому, что числа беззнаковые. Благодаря этому можно объединить 0 итерация с 1:

module divider #(parameter N = 32) ( input wire clk, input wire start, input wire [ N-1 : 0 ] divident, input wire [ N-1 : 0 ] divider, output wire [ N-1 : 0 ] quotient, output wire [ N-1 : 0 ] reminder, output wire ready ); localparam M = 2*N; reg signed [ N-1 : 0 ] r_quotient = {N{1'b0}}; assign quotient = r_quotient; reg signed [ M-1 : 0 ] divident_copy = {M{1'b0}}; reg signed [ M-1 : 0 ] divider_copy = {M{1'b0}};; wire signed [ M-1 : 0 ] w_diff = divident_copy - divider_copy; reg [ 5 : 0 ] cnt = 6'b0; assign reminder = divident_copy[N-1:0]; assign ready = cnt == 0; always@(posedge clk) if(ready && start) begin cnt <= 6'd32; r_quotient <= {N{1'b0}}; divident_copy <= {{N{1'b0}}, divident}; divider_copy <= {1'b0, divider, {N-1{1'b0}}}; end else begin cnt <= cnt - 1'b1; divider_copy <= divider_copy >> 1; if(!w_diff[63]) begin divident_copy <= w_diff; r_quotient <= {quotient[30:0], 1'b1}; end else begin r_quotient <= {quotient[30:0], 1'b0}; end end endmodule