MHB Fast Fourier Transform for Power of 3

[evinda]

Hello! (Wave)

I want to write a version of [m]FastFourierTransform(fft)[/m] for the case that $N$ is a power of $3$, seperating the input-vector into $3$ subvectors, solving the problem recursively at them and combining the solutions of the subproblems.

I have tried the following:

We assume that $N$ is a power of $3$ ($N=3^m$)

$f(x)=a_0+a_1x+ \dots +a_{N-1}x^{N-1}$

$f_0=a_0+a_3x^3+a_6x^6+ \dots +a_{N-3}x^{N-3}$

$f_1=a_1+a_4x^3+a_7x^6+ \dots +a_{N-2}x^{N-3}$

$f_2=a_2+a_5x^3+a_8x^6+ \dots +a_{N-1}x^{N-3}$

$f(x)=f_0(x)+xf_1(x)+x^2f(x)$

Now do we have to calculate the Fourier Transform of $f_0, f_1$ and $f_2$? (Thinking)Is the dicrete Fourier transform of $f$ the following? $$\sum_{m=0}^{N-1}fe^{\frac{-2\pi i km}{N}}=\sum_{m=0}^{\frac{N}{3}-1}f_{3m}e^{-\frac{2\pi ik(3m)}{N}}+\sum_{m=0}^{\frac{N}{3}}f_{3m+1}e^{-\frac{2\pi ik(3m+1)}{N}}+\sum_{m=0}^{\frac{N}{3}-1}f_{3m+2}e^{-\frac{2\pi ik(3m+2)}{N}}$$

So do we have to calculate the following sums?

• $\sum_{m=0}^{\frac{N}{3}-1}f_{3m}e^{-\frac{2\pi ik(3m)}{N}}$
  
• $\sum_{m=0}^{\frac{N}{3}}f_{3m+1}e^{-\frac{2\pi ik(3m+1)}{N}}$
  
• $\sum_{m=0}^{\frac{N}{3}-1}f_{3m+2}e^{-\frac{2\pi ik(3m+2)}{N}}$
  

 

[I like Serena]

I want to write a version of [m]FastFourierTransform(fft)[/m] for the case that $N$ is a power of $3$, seperating the input-vector into $3$ subvectors, solving the problem recursively at them and combining the solutions of the subproblems.

I have tried the following:

We assume that $N$ is a power of $3$ ($N=3^m$)

$f(x)=a_0+a_1x+ \dots +a_{N-1}x^{N-1}$

$f_0=a_0+a_3x^3+a_6x^6+ \dots +a_{N-3}x^{N-3}$

$f_1=a_1+a_4x^3+a_7x^6+ \dots +a_{N-2}x^{N-3}$

$f_2=a_2+a_5x^3+a_8x^6+ \dots +a_{N-1}x^{N-3}$

$f(x)=f_0(x)+xf_1(x)+x^2f(x)$

Now do we have to calculate the Fourier Transform of $f_0, f_1$ and $f_2$? (Thinking)

Hi! (Wink)

Yes.
Btw, I think there is no need to assume that $f$ is a polynomial.
The process works for any $f$. (Wasntme)

Is the dicrete Fourier transform of $f$ the following? $$\sum_{m=0}^{N-1}fe^{\frac{-2\pi i km}{N}}=\sum_{m=0}^{\frac{N}{3}-1}f_{3m}e^{-\frac{2\pi ik(3m)}{N}}+\sum_{m=0}^{\frac{N}{3}}f_{3m+1}e^{-\frac{2\pi ik(3m+1)}{N}}+\sum_{m=0}^{\frac{N}{3}-1}f_{3m+2}e^{-\frac{2\pi ik(3m+2)}{N}}$$

So do we have to calculate the following sums?

• $\sum_{m=0}^{\frac{N}{3}-1}f_{3m}e^{-\frac{2\pi ik(3m)}{N}}$
  
• $\sum_{m=0}^{\frac{N}{3}}f_{3m+1}e^{-\frac{2\pi ik(3m+1)}{N}}$
  
• $\sum_{m=0}^{\frac{N}{3}-1}f_{3m+2}e^{-\frac{2\pi ik(3m+2)}{N}}$
  

Yes. (Nod)

It appears your $f_{3m}$ here is different from the $f_0, f_1, f_2$ you previously defined. (Worried)
Apparently $f_{3m} = f(x_{3m}) = f(a + \frac{3m}{N}(b-a))$, where $[a,b]$ is the interval that $f$ is defined.

For the next step, I suggest to substitute $N=3M$. (Wasntme)

 

[evinda]

Hi! (Wink)

Yes.
Btw, I think there is no need to assume that $f$ is a polynomial.
The process works for any $f$. (Wasntme)

So do we assume that $f$ is any function? (Thinking)

It appears your $f_{3m}$ here is different from the $f_0, f_1, f_2$ you previously defined. (Worried)
Apparently $f_{3m} = f(x_{3m}) = f(a + \frac{3m}{N}(b-a))$, where $[a,b]$ is the interval that $f$ is defined.

For the next step, I suggest to substitute $N=3M$. (Wasntme)

Could you explain it further to me? I haven't understood it.. (Worried)

 

[evinda]

Is the general formula of the Discrete Fourier transform the following?

$$y_j= \sum_{i=0}^{n-1} a_i \omega^{ij}$$

So if $N$ is a power of $3$, do we write it as follows?$$y_j= \sum_{i=0}^{\frac{n}{3}-1} a_{3m} \omega^{(3m)j}+ \sum_{i=0}^{\frac{n}{3}-1} a_{3m+1} \omega^{(3m+1)j}+ \sum_{i=0}^{\frac{n}{3}-1} a_{3m+2} \omega^{(3m+2)j}$$

 

[evinda]

Should the algorithm look like that?
http://pastebin.com/4vuudvgt

 

[I like Serena]

So do we assume that $f$ is any function? (Thinking)

The input for an FFT is usually an array f_n = f(x_n), where $x_0, ..., x_{N-1}$ is equally spaced division of the domain of $f$.
For this, $f$ can be any function.

It's not usually a set of polynomial coefficients, as you have.
Are you supposed to use those polynomial coefficients as input instead? (Wondering)

Could you explain it further to me? I haven't understood it.. (Worried)

All I'm saying is that the discrete Fourier transform is:
$$F_k = \sum_{n=0}^{N-1} f(x_n) e^{-2\pi i kn/N}$$
where $f(x_n)$ are function values. The result $F_k$ is an array of transformed values.

Is the general formula of the Discrete Fourier transform the following?

$$y_j= \sum_{i=0}^{n-1} a_i \omega^{ij}$$

So if $N$ is a power of $3$, do we write it as follows?

$$y_j= \sum_{i=0}^{\frac{n}{3}-1} a_{3m} \omega^{(3m)j}+ \sum_{i=0}^{\frac{n}{3}-1} a_{3m+1} \omega^{(3m+1)j}+ \sum_{i=0}^{\frac{n}{3}-1} a_{3m+2} \omega^{(3m+2)j}$$

It's the right formula with a different notation, assuming $\omega = e^{-2\pi i /n}$.
But I don't think it works like that if $a_i$ are the coefficients of a polynomial.

You can rewrite it into 3 transforms though, but then the individual terms should be written as transforms.
As yet the first term contains $3m$ instead of something like $m$.
And the second term has $\omega^{(3m+1)j}$ when it should be something like $\eta^{mj}\eta$, where $\eta = e^{-2\pi i / (n/3)}$. (Wasntme)

Should the algorithm look like that?
http://pastebin.com/4vuudvgt

I believe that your formula for the transform has not been split up correctly yet. (Wasntme)

 

[evinda]

The input for an FFT is usually an array f_n = f(x_n), where $x_0, ..., x_{N-1}$ is equally spaced division of the domain of $f$.
For this, $f$ can be any function.

It's not usually a set of polynomial coefficients, as you have.
Are you supposed to use those polynomial coefficients as input instead? (Wondering)

All I'm saying is that the discrete Fourier transform is:
$$F_k = \sum_{n=0}^{N-1} f(x_n) e^{-2\pi i kn/N}$$
where $f(x_n)$ are function values. The result $F_k$ is an array of transformed values.It's the right formula with a different notation, assuming $\omega = e^{-2\pi i /n}$.
But I don't think it works like that if $a_i$ are the coefficients of a polynomial.

You can rewrite it into 3 transforms though, but then the individual terms should be written as transforms.
As yet the first term contains $3m$ instead of something like $m$.
And the second term has $\omega^{(3m+1)j}$ when it should be something like $\eta^{mj}\eta$, where $\eta = e^{-2\pi i / (n/3)}$. (Wasntme)

Is it like that?

$$F_k=\sum_{m=0}^{N-1}f(x_n)e^{-\frac{2\pi ikn}{N}} \\ =\sum_{n=0}^{\frac{N}{3}-1}f(x_{3n})e^{-\frac{2\pi ik(3n)}{N}}+\sum_{n=0}^{\frac{N}{3}-1}f(x_{3n+1})e^{-\frac{2\pi ik(3n+1)}{N}}+\sum_{n=0}^{\frac{N}{3}-1}f(x_{3n+2})e^{-\frac{2\pi ik(3n+2)}{N}} \\ =\sum_{n=0}^{\frac{N}{3}-1}f(x_{3n})e^{-\frac{2\pi ikn}{\frac{N}{3}}}+\sum_{n=0}^{\frac{N}{3}-1}f(x_{3n+1})e^{-\frac{2\pi ikn}{\frac{N}{3}}}e^{-\frac{2\pi ik}{N}}+\sum_{n=0}^{\frac{N}{3}-1}f(x_{3n+2})e^{-\frac{2\pi ikn}{\frac{N}{3}}}e^{-\frac{4\pi ik}{N}}$$

(Thinking)

I believe that your formula for the transform has not been split up correctly yet. (Wasntme)

Should it look like that: http://pastebin.com/gZKfqyYG ?

 

[I like Serena]

Is it like that?

$$F_k=\sum_{m=0}^{N-1}f(x_n)e^{-\frac{2\pi ikn}{N}} \\ =\sum_{n=0}^{\frac{N}{3}-1}f(x_{3n})e^{-\frac{2\pi ik(3n)}{N}}+\sum_{n=0}^{\frac{N}{3}-1}f(x_{3n+1})e^{-\frac{2\pi ik(3n+1)}{N}}+\sum_{n=0}^{\frac{N}{3}-1}f(x_{3n+2})e^{-\frac{2\pi ik(3n+2)}{N}} \\ =\sum_{n=0}^{\frac{N}{3}-1}f(x_{3n})e^{-\frac{2\pi ikn}{\frac{N}{3}}}+\sum_{n=0}^{\frac{N}{3}-1}f(x_{3n+1})e^{-\frac{2\pi ikn}{\frac{N}{3}}}e^{-\frac{2\pi ik}{N}}+\sum_{n=0}^{\frac{N}{3}-1}f(x_{3n+2})e^{-\frac{2\pi ikn}{\frac{N}{3}}}e^{-\frac{4\pi ik}{N}}$$

(Thinking)

Should it look like that: http://pastebin.com/gZKfqyYG ?

Yep. (Nod)

We should make a recursive call to [m]FFT3[/m] instead of [m]FFT[/m] though.
And we should give [m]N[/m] a value somewhere. (Doh)

Btw, the return value for N=2 is incorrect.
It's a good thing this cannot happen if the initial N is a power of 3.
And if it could happen, we should detect that N is not a power of 3 and return an error, since I think the algorithm won't work properly. (Nerd)

 

[evinda]

Yep. (Nod)

We should make a recursive call to [m]FFT3[/m] instead of [m]FFT[/m] though.
And we should give [m]N[/m] a value somewhere. (Doh)

Btw, the return value for N=2 is incorrect.
It's a good thing this cannot happen if the initial N is a power of 3.
And if it could happen, we should detect that N is not a power of 3 and return an error, since I think the algorithm won't work properly. (Nerd)

I changed the algorithm: http://pastebin.com/gZKfqyYG
Is it right now? (Thinking)

 

[I like Serena]

I changed the algorithm: http://pastebin.com/gZKfqyYG
Is it right now? (Thinking)

One more thing, most of the algorithm refers to $n$, but in the last line it refers to $m$, which is not consistent.
And actually, I think we might simply call them f_0, f_1, f_2, y_0, y_1, y_2, which I think is less confusing. (Wasntme)

 

[evinda]

One more thing, most of the algorithm refers to $n$, but in the last line it refers to $m$, which is not consistent.
And actually, I think we might simply call them f_0, f_1, f_2, y_0, y_1, y_2, which I think is less confusing. (Wasntme)

I changed it again: http://pastebin.com/gZKfqyYG
Is it right? (Thinking)

 

[I like Serena]

I changed it again: http://pastebin.com/gZKfqyYG
Is it right? (Thinking)

Hmm, you have:

         for i=0 to N/3-1
             y(i)=y(i)_0+y(i)_1+y(i)_2
         return y

For instance y(i)_0 looks a bit weird. It seems to me that it should be y_0(i). (Thinking)

Then again, how about making it:

         y = y_0 + y_1 + y_2
         return y

Or just:

         return y_0 + y_1 + y_2

(Wondering)

 

[evinda]

Hmm, you have:

         for i=0 to N/3-1
             y(i)=y(i)_0+y(i)_1+y(i)_2
         return y

For instance y(i)_0 looks a bit weird. It seems to me that it should be y_0(i). (Thinking)

Then again, how about making it:

         y = y_0 + y_1 + y_2
         return y

Or just:

         return y_0 + y_1 + y_2

(Wondering)

So should it look like that http://pastebin.com/gZKfqyYG ? (Thinking)

 

[I like Serena]

So should it looks like that http://pastebin.com/gZKfqyYG ? (Thinking)

That seems fine to me. (Happy)