4

On page 137 of "Understanding Digital Signal Processing" by R.G. Lyons I found that if I separate the standard DFT form:

$$X(m)=\sum_{n=0}^{N-1}x(n)e^{-j2\pi nm/N}\tag{4-11}$$

by odd and even components, with $W_{N}^{nm} = e^{-j(2\pi/N) mn}$:

$$X(m)=\sum_{n=0}^{(N/2)-1}x(2n)W_{N/2}^{nm}+W_{N}^m\sum_{n=0}^{(N/2)-1}x(2n+1)W_{N/2}^{nm}.\tag{4-15}$$ So we now have two $N/2$ summations whose results can be combined to give us the $N$-point DFT. We've reduced some of the necessary number crunching in Eq. (4-15) relative to Eq. (4-11) because the $W$ terms in the two summations of Eq. (4-15) are identical. There's a further benefit in breaking the $N$-point DFT into two parts because the upper half of the DFT outputs is easy to calculate. Consider the $X(m+N/2)$ output. If we plug $m+N/2$ in for $m$ in Eq. (4-15), then [followed by equation for $X(m + N/2)$]

I can "reduce some of the necessary number crunching, because the $W$ terms in the two summations are identical", which is not clear. It looks like changing the order of the addends does not change the sum, so the "number crunching" should be the same, isn't it?

Matt L.
  • 89,963
  • 9
  • 79
  • 179
Curious
  • 355
  • 1
  • 11

4 Answers4

8

In the original formula you need (approximately) $N^2$ complex multiplications and additions to compute all values of $X[m]$. In the second formula you have two DFTs of length $N/2$, which results in (approximately) $2(N/2)^2$ complex multiplications and additions. Furthermore, you need $N$ complex multiplications due to the term $W_N^m$ in the second sum, and you need $N$ complex additions to add the two terms. Consequently, the total number of complex multiplications and additions necessary to evaluate the second formula is (approximately) $N+2(N/2)^2=N+N^2/2$. Hence, for $N>2$ it is more efficient to split up the sum into two sums.

Note that we didn't take into account the computation of the twiddle factors $W$. We can assume that they are pre-computed and stored.

It's instructive to look at the combination of the two length $N/2$ DFT in vector notation. Let $\mathbf{X}$ denote the vector containing the elements $X[k]$, $k=0,\ldots,N-1$. Similarly, define the length $N/2$ vectors $\mathbf{X}_1$ and $\mathbf{X}_2$ containing the elements of the two length $N/2$ DFTs ($\mathbf{X}_1$ for the DFT of the even elements of $x[n]$, and $\mathbf{X}_2$ for the DFT of the odd elements of $x[n]$). The length $N$ DFT is then computed from the two length $N/2$ DFTs as follows:

$$\mathbf{X}=\left[\begin{array}{l}\mathbf{X}_1\\\mathbf{X}_1\end{array}\right]+\left[\begin{array}{l}W_N^0\\\vdots\\W_N^{N-1}\end{array}\right]\odot\left[\begin{array}{l}\mathbf{X}_2\\\mathbf{X}_2\end{array}\right]$$

where $\odot$ denotes element-wise multiplication. Note that the length $N/2$ DFTs are reused to compute the first and second halves of the length $N$ DFT.

For $N$ a power of $2$, the decimation-in-time FFT algorithm continues to split up the problem into smaller subproblems until we're left with length $2$ DFTs, which need to be combined appropriately to compute the original length $N$ DFT. The details can be found in any DSP textbook.

Matt L.
  • 89,963
  • 9
  • 79
  • 179
  • am I right: $2(N/2)^2$ complex multiplications and additions (for $W$ under the sum), $N/2$ complex multiplications (for multiplication of $W$ out of the second sum) and $N/2$ complex additions for two sums? – Curious Feb 12 '23 at 20:29
  • and another thing: here index $m$ can be defined from $0$ to $N/2$ (because $n$ is defined from $0$ to $(N/2)-1$), so the standard DFT form gives us $(N/2)^2$ multiplications, isn't it? – Curious Feb 12 '23 at 20:39
  • For the computation of the $N/2$ DFTs, the index $m$ is between $0$ and $N/2-1$. For the multiplication of the second term with $W_N^m$, the index $m$ goes up to $N-1$. That's why you get N complex additions and multiplications on top of the 2 length $N/2$ DFTs. – Matt L. Feb 13 '23 at 12:04
  • but how could it be? $(N/2)^2$ for both sums is clear, if $m$ goes from $0$ to $N/2-1$, I don't understand, why for $W_N^m$ $m$ goes from $0$ to $N-1$...this term is multiplied by the odd components, so I think it should be composed of all odd $m$ components... – Curious Feb 15 '23 at 10:52
  • 1
    @Curious: No, you're confusing indices. The odd components are the indices $n$ of the time domain signals. The index $m$ is the index of the DFT result. If you look at the derivation you should see what I mean. Also take a look at the vector equation in my answer. – Matt L. Feb 15 '23 at 12:00
  • yes, I was confusing, but now I got the issue, thank you!) by the way, why did you split $X_1$ and $X_2$ by the halves in a vector notation? – Curious Feb 15 '23 at 13:25
  • @Curious: There's no splitting. The vectors $X_1$ and $X_2$ are the length $N/2$ DFTs of the even and odd parts of the time domain sequence. So you just stack them up to get a length $N$ vector. That's the whole idea: you can reuse the length $N/2$ DFTs to compute both halves of the original length $N$ DFT. – Matt L. Feb 15 '23 at 13:37
  • I finally found the issue, with which I absolutely disagree in your answer: in your answer you consider $m$ going from $0$ to $N-1$ and NOT from $0$ to $N/2-1$, which means, that both sums give us $N^2/2$ (NOT $(N/2)^2$), so totally it will be $N^2/2+N^2/2+N = N^2 + N$, which is bigger than originally $N^2$ multiplications and summations. – Curious Feb 15 '23 at 15:01
  • 1
    @Curious: You're wrong about that. I clearly state that we got two length $N/2$ DFTs, which obviously implies that both indices in the sums go from zero to $N/2-1$. Your previous comment also shows that you're still unclear about how this works. I don't write $X=X_1+[\ldots]X_2$ because that doesn't work. The left hand side is a length $N$ vector and the first part of the right hand side is a length $N/2$ vector, whereas the product on the right hand side is invalid (length $N$ vector times length $N/2$ vector). – Matt L. Feb 15 '23 at 15:08
  • ok, as I understand the idea now: we just calculate the coefficients for one half of the FT ($(N/2)-1$) and then change the sign for the second sum; but in this case your vector notation looks wrong for me, because for the second part of the spectrum ($X(m+N/2)$) the second sum has the negative sign (before $W_N^m$)... – Curious Feb 15 '23 at 15:18
  • 1
    The sign is taken into account by the factors $W_N^m$ because in that sum (and only there) does the index $m$ go from $0$ to $N-1$. You could let $m$ go from $0$ to $N/2-1$ and change sign for the second half, but I chose to write it that way. Anyway, the sum is written correctly. I hope you agree. – Matt L. Feb 15 '23 at 15:27
  • 1
    now (I hope) the idea looks absolutely clear for me) – Curious Feb 15 '23 at 15:45
  • I just was thinking about this issue again, what if we'll calculate the spectrum components $m$ only from $0$ to $N/2-1$ (because we clearly know, that we'll get the "mirrored" samples relative to $N/2$); in this case the standard DFT requires $N^2/2$ complex multiplications and summations and it looks like the use of FFT makes the calculation even worse (because it requires $N^2/2+N$ complex multiplications and summations)... – Curious Mar 04 '23 at 15:43
  • 1
    @Curious: First, that symmetry holds only for real-valued sequences, and the DFT can generally be applied to complex-valued sequences. Second, the complexity of the FFT is not what you wrote. That step discussed in your question and my answer is just a single step of the FFT algorithm. Read the last paragraph of my answer. The complexity of the FFT is proportional to $N\log N$. – Matt L. Mar 04 '23 at 20:38
5

Your computation is reduced because $N$ is even, $(N/2)$ is an integer and $W_N^{(N/2)}=-1$, $W_N^2=W_{N/2}$, which implies $W_{(N/2)}^{(N/2)}=1$, so that the partial sums can be seen as DFT of half the length. This allows to decompose the DFT sums of length $M$ into parts that are DFT sums of length $N/2$, \begin{alignat}{1} X(m+(N/2))&=&&\sum_{n=0}^{(N/2)-1}x(2n)W_{N/2}^{nm}\\&&-W_{N}^m&\sum_{n=0}^{(N/2)-1}x(2n+1)W_{N/2}^{nm} \end{alignat} This computes an extra value using the same partial sums as used for $X(m)$.

If this is recursively continued, one gets the complexity equation $$T(N)=2T(N/2)+c\cdot(N/2),$$ where $T(N)$ is the time or complexity of transforming an array of size $N$ and $c$ some small constant collecting the operations outside the sums. This gives $$T(N)=(c/2)⋅N\log_2(N)$$ if $N$ is a power of $2$. There are other FFT algorithms for small factors other than 2 that allow to extend this pattern to general sizes $N$.

Lutz Lehmann
  • 1,185
  • 5
  • 9
  • thank you for your answer, but I'm not familiar with complexity equations and don't actually know, what do $T(N)$ and $c(N/2)$ mean(( – Curious Feb 15 '23 at 14:01
3

If you calculate the DFT naively, you do a lot of duplicate work. Here is the DFT formula again

$$X(m)=\sum_{n=0}^{N-1}x(n)\omega^{nm}$$ where ω is a Nth complex root of unity, meaning $$\begin{align*} \omega=e^\frac{-j2\pi }{N} && \text{with} && \omega^N = 1 && \text{and} && \omega^{nm} = \omega^{N+nm} = \omega^{2N+nm} = \space ...\\ \end{align*}$$

That means you have lots of {n,m} tuples that result in the same power of ω. For N=20, for example, there is:

$$ \omega^{3\cdot5} = \omega^{5\cdot3} = \omega^{5\cdot7} = \omega^{7\cdot5} = \omega^{5\cdot11} = \omega^{11\cdot5}= \omega^{3\cdot25} =\space...\space= \omega^{kN+15} = \omega^{15} $$

All FFT algorithms re-use intermediate results to avoid duplication. There's a lot of duplication (N² {n,m} tuples, but only N distinct powers of ω) and that's why the non-naive algorithms are much faster.

Rainer P.
  • 531
  • 2
  • 6
  • cool! I wish I could accept several answers(( – Curious Feb 15 '23 at 15:43
  • by the way, what is the need of calculating the samples above $N/2$? if we clearly know, that another half of the samples will be mirrored relative to $N/2$ sample? – Curious Mar 04 '23 at 16:02
  • If that were the main reason, you could reduce the complexity by using an array for the powers of $\omega$. – Lutz Lehmann Dec 07 '23 at 10:42
  • @LutzLehmann - An array for the powers of ω is not enough. The fact that you encounter the same powers of ω multiple times also means that you encounter the same partial sums multiple times. You want to reuse partial sums in a hierarchical manner to avoid duplicate work. – Rainer P. Dec 07 '23 at 13:17
  • Yes, see my answer that I touched up alongside the comment. It is indeed the number of multiplications of unit root powers and data values that one needs to reduce, not only avoid the redundancy in the unit-root-power computations. – Lutz Lehmann Dec 07 '23 at 13:28
2

Looking just at the formulas and not jumping ahead deeper into FFT, I got to agree with you that there are no direct savings in computation due to the coefficients being identical in the two sums. Let's try to combine two multiplications by identical coefficients into a single multiplication:

$$\begin{eqnarray}X(m)&=&\sum_{n=0}^{(N/2)-1}x(2n)W_{N/2}^{nm}+W_{N}^m\sum_{n=0}^{(N/2)-1}x(2n+1)W_{N/2}^{nm}\\ &=&\sum_{n=0}^{(N/2)-1}x(2n)W_{N/2}^{nm}+\sum_{n=0}^{(N/2)-1}W_{N}^mx(2n+1)W_{N/2}^{nm}\\ &=&\sum_{n=0}^{(N/2)-1}\left(x(2n)W_{N/2}^{nm}+W_{N}^mx(2n+1)W_{N /2}^{nm}\right)\\ &=&\sum_{n=0}^{(N/2)-1}\left(x(2n)+W_{N}^mx(2n+1)\right)W_{N /2}^{nm},\end{eqnarray}$$

which for a single value of $m$ would be going from $N/2 + N/2 + 1 = N + 1$ multiplications and $(N/2 - 1) + (N/2 - 1) + 1 = N-1$ additions on the first row to $2\,N/2 = N$ multiplications and $N/2 + N/2 - 1 = N-1$ additions on the last row. But that's the same as $N$ multiplications and $N-1$ additions in the original formula that we were comparing against:

$$X(m)=\sum_{n=0}^{N-1}x(n)W_N^{nm}.$$

I don't know if Rick just refers to needing to (pre)calculate less values of $W$.

Olli Niemitalo
  • 13,491
  • 1
  • 33
  • 61
  • 1
    Maybe I'm missing something, but did you account for the fact that the index $m$ in the two sums of your first expression is only evaluated from $m=0$ to $m=N/2-1$? So you have two length $N/2$ DFTs with (approx.) $N^2/4$ complex additions and multiplications each, plus the overhead for adding the two terms and multiplying the second term with $W_N^m$ (this time for all $m$ up to $N-1$). This should result in the saving I mentioned in my answer. – Matt L. Feb 13 '23 at 12:08
  • @MattL. For me $m$ still goes up to $N-1$, because that's how it still is at the point in the book that the question was about. I have added a quote from the book to the question to clarify. – Olli Niemitalo Feb 13 '23 at 13:01
  • But he actually continues to show that it is sufficient to let $m$ go from $0$ to $N/2-1$ in the two sums. The proof could be more straightforward, but on page 138 he finally arrives at that conclusion. "Of course, $m$ goes from $0$ to $N/2-1$ in Eq. (4-20) ...". – Matt L. Feb 13 '23 at 14:11
  • 1
    @MattL. I'm just anal about separating 1) a purported reduction in necessary number crunching due to the $W$ terms being equal in the two summations of Eq. (5-14) and 2) "a further benefit in breaking the $N$-point DFT into two parts" which is where the computational savings come from and which you and Lutz Lehmann have elaborated on. – Olli Niemitalo Feb 13 '23 at 14:36
  • 1
    @OlliNiemitalo, you're the one who got the issue of my question!) I don't think that precalculating of $W$ will make the difference compared to the standard DFT form (because it is also depend on both $m$ and $n$). – Curious Feb 13 '23 at 15:10