Покажите, как сделать FFT вручную

27

Скажем, у вас есть два полинома: и .3+x2x2+2

Я пытаюсь понять, как БПФ помогает нам умножить эти два полинома. Однако я не могу найти какие-либо разработанные примеры. Может кто-нибудь показать мне, как алгоритм FFT умножит эти два полинома. (Примечание: в этих многочленах нет ничего особенного, но я хотел, чтобы это было просто, чтобы было легче следовать.)

Я посмотрел на алгоритмы в псевдокоде, но все они, похоже, имеют проблемы (не указывайте, какой должен быть ввод, неопределенные переменные). И, что удивительно, я не могу найти, где кто-нибудь на самом деле прошел (вручную) пример умножения многочленов с использованием БПФ.

Lars
источник
2
Википедия хранит это прекрасное изображение для целочисленного умножения через FFT, но я думаю, что еще более явный пошаговый шаг может быть полезным.
Реал Слав

Ответы:

27

Предположим, что мы используем четвертые корни единицы, что соответствует замене 1,i,1,i на x . Мы также используем прореживание по времени, а не прореживание по частоте в алгоритме FFT. (Мы также без проблем применяем операцию обращения битов.)

Чтобы вычислить преобразование первого полинома, мы начнем с записи коэффициентов:

3,1,0,0.
Преобразование Фурье четных коэффициентов 3,0 равно 3,3 , а нечетных коэффициентов 1,0 равно 1,1 . (Это преобразование просто a,ba+b,ab .) Следовательно, преобразование первого полинома равно
4,3+i,2,3i.
Это получается с использованиемX0,2=E0±O0 ,X1,3=E1iO1 . (Из расчета тиддл-фактора).

Давайте сделаем то же самое для второго полинома. Коэффициенты равны

2,0,2,0.
Четные коэффициенты 2,2 преобразуются в 4,0 , а нечетные коэффициенты 0,0 преобразуются в 0,0 . Следовательно, преобразование второго полинома равно
4,0,4,0.

Мы получаем преобразование Фурье полинома произведений, умножая два преобразования Фурье поточечно:

16,0,8,0.
Осталось вычислить обратное преобразование Фурье. Четные коэффициенты 16,8 обратного преобразования в 12,4 и нечетные коэффициенты 0,0 обратного преобразования в 0,0 . (Обратное преобразование: x,y(x+y)/2,(xy)/2 ) Следовательно, преобразование полинома произведений равно
6,2,6,2.
Это получается с помощьюX0,2=(E0±O0)/2 ,X1,3=(E1iO1)/2 . Мы получили желаемый ответ
(3+x)(2+2x2)=6+2x+6x2+2x3.

Юваль Фильмус
источник
как ты добрался до 6,2 6,2?
Ларс
Я дал формулы: , X 1 , 3 = ( E 1i O 1 ) / 2 , где E 0 , E 1 ( O 1 , O 2 ) обратное преобразование четных (нечетных) коэффициентов, полученное по формуле x , y ( x + y )X0,2=(E0±O2)/2X1,3=(E1iO1)/2E0,E1O1,O2 . Пожалуйста, посмотрите на ответ еще раз - все расчеты есть. x,y(x+y)/2,(xy)/2
Юваль Фильмус
Почему вы используете четные коэффициенты дважды? 3,3 -> 3,3,3,3. -> 3 + 1, 3-я, 3 + -1,3 - я?
Aage Torleif
Как эти формулы для и X 1 , 3 распространяются на более высокие степени? Знаки плюс / минус просто продолжают переворачиваться? Например, что это будет за X 0 , 2 , 4 ? X0,2X1,3X0,2,4
Бобби Ли,
@BobbyLee Я призываю вас прочитать литературу по FFT.
Юваль Фильмус
7

Определите полиномы, где deg(A) = qи deg(B) = p. deg(C) = q + p.

В этом случае deg(C) = 1 + 2 = 3.

A=3+xB=2x2+2C=AB=?

Мы легко можем найти C за время O(n2) путем умножения коэффициентов методом грубой силы. Применяя БПФ (и обратное БПФ), мы могли бы достичь этого за O(nlog(n)) времени. Явное:

  1. Преобразуйте представление коэффициента A и B в его представление значения. Этот процесс называется оценкой . Выполнение «Разделяй и властвуй» (D & C) для этого займет время O(nlog(n)) .
  2. Умножьте многочлены по компонентам в их представлении значения. Это возвращает представление значения C = A * B. Это займет O(n) время.
  3. Инвертируйте C, используя обратное БПФ, чтобы получить C в его представлении коэффициента. Этот процесс называется интерполяцией и также занимает время O(nlog(n)) .

Continuing along, we represent each polynomial as a vector whose value are its coefficients. We pad the vector with 0's up to the smallest power of two, n=2k,ndeg(C). Thus n=4. Choosing a power of two provides us a way to recursively apply our divide-and-conquer algorithm.

A=3+x+0x2+0x3a=[3,1,0,0]B=2+0x+2x+0x3b=[2,0,2,0]

Let A,B be the value representation of A and B, respectively. Notice that FFT (Fast Fourier Transform) is a linear transformation (linear map) and can be represented as a matrix, M. Thus

A=MaB=Mb

We define M=Mn(ω) where ω is complex roots nth complex roots of unity. Notice n = 4, in this example. Also notice that the entry in the jth row and kth column is ωnjk . See more about the DFT matrix here

M4(w)=[111...11ω1ω2...ωn11ω2ω4...............ωjk...1ωn1ω2(n1)...ω(n1)(n1)]=[11111ωω2ω31ω2ω4ω61ω3ω6ω9]

Given the ω4=4th roots of unity, we have the ordered set equality:

{ω0,ω1,ω2,ω3,ω4,ω5,...}={1,i,1,i,1,i,...}

This can be visualized as iterating thru roots of the unit circle in the counter-clockwise direction.

Also, notice the mod n nature, i.e. ω6=ω6modn=ω2=1 and i=ω3=ω3+n

To complete step 1 (evaluation) we find A,B by performing

A=Ma=[11111ωω2ω31ω2ω4ω61ω3ω6ω9][3100]=[3+13+1ω3+ω23+ω3]=[43+i23i]B=Mb=[11111ωω2ω31ω2ω4ω61ω3ω6ω9][2020]=[2+22+2ω22+2ω42+2ω6]=[4040]

This step can be achieved using D&C algorithms (beyond the scope of this answer).

Multiply AB component-wise (step 2)

AB=[43+i23i][4040]=[16080]=C

Finally, the last step is to represent C' into coefficients. Notice

C=McM1C=M1Mcc=M1C

Notice Mn1=1nMn(ω1)1 and ωj=ωn/2+j.

Mn1=14[11111ω1ω2ω31ω2ω4ω61ω3ω6ω9]=14[11111i1i11111i1i]

ωj can be visualized as iterating thru roots of the unit circle in the clockwise direction.

{ω0,ω1,ω2,ω3,ω4,ω5,...}={1,i,1,i,1,i,...}

Also, it is true that, given the nth root of unity, the equality ωj=ωnj holds. (Do you see why?)

Then,

c=M1C=1nMn(w1)=14[11111i1i11111i1i][16080]=[(16+8)/4(168)/4(16+8)/4(168)/4]=[6262]

Thus, we get the polynomial

C=AB=6+2x+6x2+2x3
1: Inversion Formula pg 73, Algorithms by Dasgupta et. al. (C) 2006

j__gt
источник