View source for Discrete Fourier transform over a ring

{{Short description|Generalisation of Fourier transform to any ring}}
{{Fourier transforms}}

In [[mathematics]], the '''discrete Fourier transform over a ring''' generalizes the [[discrete Fourier transform]] (DFT), of a function whose values are commonly [[complex number]]s, over an arbitrary [[ring (mathematics)|ring]].

==Definition==
Let {{mvar|R}} be any [[ring (mathematics)|ring]], let <math>n\geq 1</math> be an integer, and let <math>\alpha \in R</math> be a [[Principal root of unity|principal]] ''n''th root of unity, defined by:<ref name="furer">Martin Fürer, "[https://web.archive.org/web/20130425232048/http://www.cse.psu.edu/~furer/Papers/mult.pdf Faster Integer Multiplication]", STOC 2007 Proceedings, pp. 57&ndash;66. Section 2: The Discrete Fourier Transform.</ref>

{{NumBlk|:| <math>
\begin{align}
& \alpha^n = 1 \\
& \sum_{j=0}^{n-1} \alpha^{jk} = 0 \text{ for } 1 \leq k < n 
\end{align}
</math>|{{EquationRef|1}}}}

The discrete Fourier transform maps an [[n-tuple|''n''-tuple]] <math>(v_0,\ldots,v_{n-1})</math> of elements of {{mvar|R}} to another ''n''-tuple <math>(f_0,\ldots,f_{n-1})</math> of elements of {{mvar|R}} according to the following formula:

{{NumBlk|:|<math>f_k = \sum_{j=0}^{n-1} v_j\alpha^{jk}.</math>|{{EquationRef|2}}}}

By convention, the tuple <math>(v_0,\ldots,v_{n-1})</math> is said to be in the ''time domain'' and the index {{mvar|j}} is called ''time''. The tuple <math>(f_0,\ldots,f_{n-1})</math> is said to be in the ''frequency domain'' and the index {{mvar|k}} is called ''frequency''. The tuple <math>(f_0,\ldots,f_{n-1})</math> is also called the ''[[frequency spectrum|spectrum]]'' of <math>(v_0,\ldots,v_{n-1})</math>. This terminology derives from the applications of Fourier transforms in [[signal processing]].

If {{mvar|R}} is an [[integral domain]] (which includes [[field (mathematics)|fields]]), it is sufficient to choose <math>\alpha</math> as a [[primitive root of unity|primitive ''n''th root of unity]], which replaces the condition ({{EquationNote|1}}) by:<ref name="furer"/>

:<math>\alpha^{k} \ne 1</math> for <math>1 \leq k < n</math>

{{Math proof|1=Take <math>\beta = \alpha^k</math> with <math>1 \leq k < n</math>. Since <math>\alpha^n=1</math>, <math>\beta^n=(\alpha^n)^k=1</math>, giving:
:<math>\beta^n-1 = (\beta-1)\left(\sum_{j=0}^{n-1} \beta^j\right) = 0</math>
where the sum matches ({{EquationNote|1}}). Since <math>\alpha</math> is a primitive root of unity, <math>\beta - 1 \ne 0</math>. Since {{mvar|R}} is an integral domain, the sum must be zero. ∎}}

Another simple condition applies in the case where ''n'' is a power of two: ({{EquationNote|1}}) may be replaced by <math>\alpha^{n/2} = -1</math>.<ref name="furer"/>

==Inverse==

The inverse of the discrete Fourier transform is given as:

{{NumBlk|:|<math>v_j = \frac{1}{n}\sum_{k=0}^{n-1} f_k\alpha^{-jk}.</math>|{{EquationRef|3}}}}

where <math>1/n</math> is the multiplicative inverse of {{mvar|n}} in {{mvar|R}} (if this inverse does not exist, the DFT cannot be inverted).

{{Math proof|1=Substituting ({{EquationNote|2}}) into the right-hand-side of ({{EquationNote|3}}), we get

:<math>
\begin{align}
& \frac{1}{n}\sum_{k=0}^{n-1} f_k\alpha^{-jk} \\
= {} & \frac{1}{n}\sum_{k=0}^{n-1}\sum_{j'=0}^{n-1} v_{j'}\alpha^{j'k}\alpha^{-jk} \\
= {} & \frac{1}{n}\sum_{j'=0}^{n-1} v_{j'} \sum_{k=0}^{n-1}\alpha^{(j'-j)k}.
\end{align}
</math>

This is exactly equal to <math>v_j</math>, because 
<math>\sum_{k=0}^{n-1}\alpha^{(j'-j)k}=0</math> when <math>j'\neq
j</math> (by ({{EquationNote|1}}) with <math>k=j'-j</math>), and 
<math>\sum_{k=0}^{n-1}\alpha^{(j'-j)k}=n</math> when <math>j'=j</math>. ∎}}

==Matrix formulation==

Since the discrete Fourier transform is a [[linear operator]], it can be described by [[matrix multiplication]]. In matrix notation, the discrete Fourier transform is expressed as follows:
:<math>
\begin{bmatrix}f_0\\f_1\\\vdots\\f_{n-1}\end{bmatrix}
= \begin{bmatrix}
1&1&1&\cdots &1 \\
1&\alpha&\alpha^2&\cdots&\alpha^{n-1} \\
1&\alpha^2&\alpha^4&\cdots&\alpha^{2(n-1)}\\
\vdots&\vdots&\vdots&\ddots&\vdots\\
1&\alpha^{n-1}&\alpha^{2(n-1)}&\cdots&\alpha^{(n-1)(n-1)}\\
\end{bmatrix}
\begin{bmatrix}v_0\\v_1\\\vdots\\v_{n-1}\end{bmatrix}.
</math>
The matrix for this transformation is called the [[DFT matrix]].

Similarly, the matrix notation for the inverse Fourier transform is 
:<math>
\begin{bmatrix}v_0\\v_1\\\vdots\\v_{n-1}\end{bmatrix}
= \frac{1}{n}\begin{bmatrix}
1&1&1&\cdots &1 \\
1&\alpha^{-1}&\alpha^{-2}&\cdots&\alpha^{-(n-1)} \\
1&\alpha^{-2}&\alpha^{-4}&\cdots&\alpha^{-2(n-1)}\\
\vdots&\vdots&\vdots&\ddots&\vdots\\
1&\alpha^{-(n-1)}&\alpha^{-2(n-1)}&\cdots&\alpha^{-(n-1)(n-1)}
\end{bmatrix}
\begin{bmatrix}f_0\\f_1\\\vdots\\f_{n-1}\end{bmatrix}.
</math>

==Polynomial formulation==

Sometimes it is convenient to identify an {{mvar|n}}-tuple <math>(v_0,\ldots,v_{n-1})</math> with a formal polynomial

:<math>p_v(x) = v_0 + v_1x + v_2x^2 + \cdots + v_{n-1}x^{n-1}. \, </math>

By writing out the summation in the definition of the discrete Fourier transform ({{EquationNote|2}}), we obtain:

:<math>f_k = v_0 + v_1\alpha^{k} + v_2\alpha^{2k}  + \cdots + v_{n-1}\alpha^{(n-1)k}. \, </math>

This means that <math>f_k</math> is just the value of the polynomial <math>p_v(x)</math> for <math>x=\alpha^k</math>, i.e.,

{{NumBlk|:|<math>f_k = p_v(\alpha^k).\,</math>|{{EquationRef|4}}}}

The Fourier transform can therefore be seen to relate the ''coefficients'' and the ''values'' of a polynomial: the coefficients are in the time-domain, and the values are in the frequency domain. Here, of course, it is important that the polynomial is evaluated at the {{mvar|n}}th roots of unity, which are exactly the powers of <math>\alpha</math>.

Similarly, the definition of the inverse Fourier transform ({{EquationNote|3}}) can be written:

{{NumBlk|:|<math>v_j = \frac{1}{n}(f_0 + f_1\alpha^{-j} + f_2\alpha^{-2j}  + \cdots + f_{n-1}\alpha^{-(n-1)j}).</math>|{{EquationRef|5}}}}

With

:<math>p_f(x) = f_0 + f_1x + f_2x^2 + \cdots + f_{n-1}x^{n-1},</math>

this means that

:<math>v_j = \frac{1}{n}p_f(\alpha^{-j}).</math>

We can summarize this as follows: if the ''values'' of <math>p_v(x)</math> are the ''coefficients'' of <math>p_f(x)</math>, then the ''values'' of <math>p_f(x)</math> are the ''coefficients'' of <math>p_v(x)</math>, up to a scalar factor and reordering.<ref>R. Lidl and G. Pilz. Applied Abstract Algebra, 2nd edition. Wiley, 1999, pp. 217–219.</ref>

==Special cases==

===Complex numbers===

If <math>F={\mathbb C}</math> is the field of complex numbers, then the <math>n</math>th roots of unity can be visualized as points on the [[unit circle]] of the [[complex plane]]. In this case, one usually takes

:<math>\alpha=e^{\frac{-2\pi i}{n}},</math>

which yields the usual formula for the [[discrete Fourier transform|complex discrete Fourier transform]]:

:<math>f_k = \sum_{j=0}^{n-1} v_j e^{\frac{-2\pi i}{n}jk}.</math>

Over the complex numbers, it is often customary to normalize the formulas for the DFT and inverse DFT by using the scalar factor <math>\frac{1}{\sqrt{n}}</math> in both formulas, rather than <math>1</math> in the formula for the DFT and <math>\frac{1}{n}</math> in the formula for the inverse DFT. With this normalization, the DFT matrix is then unitary. 
Note that <math>\sqrt{n}</math> does not make sense in an arbitrary field.

===Finite fields===
If <math>F=\mathrm{GF}(q)</math> is a [[finite field]], where {{mvar|q}} is a [[prime number|prime]] power, then the existence of a primitive {{mvar|n}}th root automatically implies that {{mvar|n}} [[divides]] <math>q-1</math>, because the [[multiplicative order]] of each element must divide the size of the [[multiplicative group]] of {{mvar|F}}, which is <math>q-1</math>. This in particular ensures that <math>n=\underbrace{1+1+\cdots+1}_{n\ \rm times}</math> is invertible, so that the notation <math>\frac{1}{n}</math> in ({{EquationNote|3}}) makes sense.

An application of the discrete Fourier transform over <math>\mathrm{GF}(q)</math> is the reduction of [[Reed–Solomon code]]s to [[BCH code]]s in [[coding theory]]. Such transform can be carried out efficiently with proper fast algorithms, for example, [[cyclotomic fast Fourier transform]].

==== Polynomial formulation without n<sup>th</sup> root ====

Suppose <math>F=\mathrm{GF}(p)</math>. If <math>p\nmid n</math>, it may be the case that <math>n \nmid p-1</math>. This means we cannot find an <math>n^{th}</math> root of unity in <math>F</math>. We may view the Fourier transform as an isomorphism <math>\mathrm{F}[C_n]=\mathrm{F}[x]/(x^n-1) \cong \bigoplus_i \mathrm{F}[x]/(P_i(x))</math> for some polynomials <math>P_i(x)</math>, in accordance with [[Maschke's theorem]]. The map is given by the [[Chinese remainder theorem]], and the inverse is given by applying [[Bézout's identity]] for polynomials.<ref>{{cite web | url=https://github.com/jacksonwalters/dft-finite-groups | title=Jacksonwalters/DFT-finite-groups | website=[[GitHub]] }}</ref>

<math>x^n-1=\prod_{d|n} \Phi_d(x)</math>, a product of cyclotomic polynomials. Factoring <math>\Phi_d(x)</math> in <math>F[x]</math> is equivalent to factoring the prime ideal <math>(p)</math> in <math>\mathrm{Z}[\zeta] = \mathrm{Z}[x]/(\Phi_d(x))</math>. We obtain <math>g</math> polynomials <math>P_1 \ldots P_g</math> of degree <math>f</math> where <math>fg = \varphi(d)</math> and <math>f</math> is the order of <math>p \text{ mod } d</math>.

As above, we may extend the base field to <math>\mathrm{GF}(q)</math> in order to find a primitive root, i.e. a splitting field for <math>x^n-1</math>. Now <math>x^n-1 = \prod_k (x-\alpha^{k})</math>, so an element <math>\sum_{j=0}^{n-1} v_j x^j \in F[x]/(x^n-1)</math> maps to <math>\sum_{j=0}^{n-1} v_j x^j \mod (x-\alpha^k) \equiv \sum_{j=0}^{n-1} v_j (\alpha^k)^j</math> for each <math>k</math>.

==== When p divides n ====

When <math>p|n</math>, we may still define an <math>F_p</math>-linear isomorphism as above. Note that <math>(x^n-1)=(x^m-1)^{p^s}</math> where <math>n=m p^s</math> and <math>p \nmid m</math>. We apply the above factorization to <math>x^m-1</math>, and now obtain the decomposition <math>F[x]/(x^n-1) \cong \bigoplus_i F[x]/(P_i(x)^{p^s})</math>. The modules occurring are now indecomposable rather than irreducible.

==== Order of the DFT matrix ====

Suppose <math>p \nmid n</math> so we have an <math>n^{th}</math> root of unity <math>\alpha</math>. Let <math>A</math> be the above DFT matrix, a Vandermonde matrix with entries <math>A_{ij}=\alpha^{ij}</math> for <math>0 \le i,j < n</math>. Recall that <math>\sum_{j=0}^{n-1} \alpha^{(k-l)j} = n\delta_{k,l}</math> since if <math>k=l</math>, then every entry is 1. If <math>k \ne l</math>, then we have a geometric series with common ratio <math>\alpha^{k-l}</math>, so we obtain <math>\frac{1-\alpha^{n(k-l)}}{1-\alpha^{k-l}}</math>. Since <math>\alpha^n=1</math> the numerator is zero, but <math>k-l \ne 0</math> so the denominator is nonzero.

First computing the square, <math>(A^2)_{ik} = \sum_{j=0}^{n-1} \alpha^{j(i+k)} = n\delta_{i,-k}</math>. Computing <math>A^4=(A^2)^2</math> similarly and simplifying the deltas, we obtain <math>(A^4)_{ik}=n^2\delta_{i,k}</math>. Thus, <math>A^4=n^2 I_n</math> and the order is <math>4\cdot \text{ord}(n^2)</math>.

==== Normalizing the DFT matrix ====

In order to align with the complex case and ensure the matrix is order 4 exactly, we can normalize the above DFT matrix <math>A</math> with <math>\frac{1}{\sqrt{n}}</math>. Note that though <math>\sqrt{n}</math> may not exist in the splitting field <math>F_q</math> of <math>x^n-1</math>, we may form a quadratic extension <math>F_{q^2} \cong F_q[x]/(x^2-n)</math> in which the square root exists. We may then set <math>U=\frac{1}{\sqrt{n}}A</math>, and <math>U^4=I_n</math>.

==== Unitarity ====

Suppose <math>p \nmid n</math>. One can ask whether the DFT matrix is [[unitary group#finite fields|unitary over a finite field]]. If the matrix entries are over <math>F_q</math>, then one must ensure <math>q</math> is a perfect square or extend to <math>F_{q^2}</math> in order to define the order two automorphism <math>x \mapsto x^q</math>. Consider the above DFT matrix <math>A_{ij}=\alpha^{ij}</math>. Note that <math>A</math> is symmetric. Conjugating and transposing, we obtain <math>A_{ij}^{*}=\alpha^{qji}</math>.

<math>(AA^*)_{ik} = \sum_{j=0}^{n-1}\alpha^{j(i+qk)} = n\delta_{i,-qk}</math>

by a similar geometric series argument as above. We may remove the <math>n</math> by normalizing so that <math>U = \frac{1}{\sqrt{n}}A</math> and <math>(UU^*)_{ik} = \delta_{i,-qk}</math>. Thus <math>U</math> is unitary iff <math>q \equiv -1 \,(\text{mod} \, n)</math>. Recall that since we have an <math>n^{th}</math> root of unity, <math>n|q^2-1</math>. This means that <math>q^2 - 1 \equiv (q+1)(q-1) \equiv 0 \,(\text{mod} \, n)</math>. Note if <math>q</math> was not a perfect square to begin with, then <math>n|q-1</math> and so <math>q \equiv 1 \, (\text{mod} \, n)</math>.

For example, when <math>p=3, n=5</math> we need to extend to <math>q^2=3^4</math> to get a 5th root of unity. <math>q=9 \equiv -1 \, (\text{mod} \, 5)</math>.

For a nonexample, when <math>p=3, n=8</math> we extend to <math>F_{3^2}</math> to get an 8th root of unity. <math>q^2=9</math>, so <math>q \equiv 3 \,(\text{mod} \, 8)</math>, and in this case <math>q+1 \not\equiv 0</math> and <math>q-1 \not\equiv 0</math>. <math>UU^*</math> is a square root of the identity, so <math>U</math> is not unitary.

==== Eigenvalues of the DFT matrix ====

When <math>p \nmid n</math>, we have an <math>n^{th}</math> root of unity <math>\alpha</math> in the splitting field <math>F_q \cong F_p[x]/(x^n-1)</math>. Note that the characteristic polynomial of the above DFT matrix may not split over <math>F_q</math>. 
The DFT matrix is order 4. We may need to go to a further extension <math>F_{q'}</math>, the splitting extension of the characteristic polynomial of the DFT matrix, which at least contains fourth roots of unity. If <math>a</math> is a generator of the multiplicative group of <math>F_{q'}</math>, then the eigenvalues are <math>\{\pm 1,\pm a^{(q'-1)/4}\}</math>, in exact analogy with the complex case. They occur with some nonnegative multiplicity.

===Number-theoretic transform===
The '''number-theoretic transform (NTT)'''<ref>{{Cite journal|last1=Agarwal|first1=R.|last2=Burrus|first2=C.|date=April 1974|title=Fast Convolution using fermat number transforms with applications to digital filtering|url=https://ieeexplore.ieee.org/document/1162555|journal=IEEE Transactions on Acoustics, Speech, and Signal Processing|language=en|volume=22|issue=2|pages=87–97|doi=10.1109/TASSP.1974.1162555|issn=0096-3518}}</ref> is obtained by specializing the discrete Fourier transform to <math>F={\mathbb Z}/p</math>, the [[modular arithmetic|integers modulo a prime {{mvar|p}}]]. This is a [[finite field]], and primitive {{mvar|n}}th roots of unity exist whenever {{mvar|n}} divides <math>p-1</math>, so we have <math>p=\xi n+1</math> for a positive integer {{mvar|&xi;}}. Specifically, let <math>\omega</math> be a primitive <math>(p-1)</math>th root of unity, then an {{mvar|n}}th root of unity <math>\alpha</math> can be found by letting <math>\alpha=\omega^{\xi}</math>.

e.g. for <math>p=5</math>, <math>\alpha = 2</math>

:<math>\begin{align}2^{1}&=2 \pmod 5\\2^{2}&=4 \pmod 5\\2^{3}&=3 \pmod 5\\2^{4}&=1 \pmod 5\end{align}</math>

when <math>N=4</math>
:<math>
\begin{bmatrix}
F(0) \\
F(1) \\
F(2) \\
F(3) \end{bmatrix}
=
 \begin{bmatrix}
1 & 1 & 1 & 1 \\
1 & 2 & 4 & 3 \\
1 & 4 & 1 & 4 \\
1 & 3 & 4 & 2 \end{bmatrix}
\begin{bmatrix}
f(0) \\
f(1) \\
f(2) \\
f(3) \end{bmatrix}
</math>
 
The number theoretic transform may be meaningful in the [[ring (mathematics)|ring]] <math>\mathbb{Z}/m</math>, even when the modulus {{mvar|m}} is  not prime, provided a principal root of order {{mvar|n}} exists. Special cases of the number theoretic transform such as the Fermat Number Transform ({{math|1=''m'' = 2<sup>''k''</sup>+1}}), used by the [[Schönhage–Strassen algorithm]], or Mersenne Number Transform<ref>{{Cite journal|last=Rader|first=C.M.|date=December 1972|title=Discrete Convolutions via Mersenne Transforms|url=https://ieeexplore.ieee.org/document/1672090|journal=IEEE Transactions on Computers|volume=C-21|issue=12|pages=1269–1273|doi=10.1109/T-C.1972.223497|s2cid=1939809 |issn=0018-9340}}</ref> ({{math|1=''m'' = 2<sup>''k''</sup>&nbsp;&minus;&nbsp;1}}) use a composite modulus.

===Discrete weighted transform===
The '''discrete weighted transform (DWT)''' is a variation on the discrete Fourier transform over arbitrary rings involving [[weight function|weighting]] the input before transforming it by multiplying elementwise by a weight vector, then weighting the result by another vector.<ref>{{Citation | last1 = Crandall | first1 = Richard | last2 = Fagin | first2 = Barry | title = Discrete weighted transforms and large-integer arithmetic | journal = Mathematics of Computation | volume = 62 | issue = 205 | year = 1994 | pages = 305&ndash;324 | url = http://www.faginfamily.net/barry/Papers/Discrete%20Weighted%20Transforms.pdf | doi=10.2307/2153411| jstor = 2153411 | doi-access = free }}</ref> The [[Irrational base discrete weighted transform]] is a special case of this.

==Properties==

Most of the important attributes of the [[discrete Fourier transform|complex DFT]], including the inverse transform, the [[convolution theorem]], and most [[fast Fourier transform]] (FFT) algorithms, depend only on the property that the kernel of the transform is a principal root of unity. These properties also hold, with identical proofs, over arbitrary rings. In the case of fields, this analogy can be formalized by the [[field with one element]], considering any field with a primitive ''n''th root of unity as an algebra over the extension field <math>\mathbf{F}_{1^n}.</math>{{Clarify|reason=I don't understand|date=November 2018}}

In particular, the applicability of <math>O(n \log n)</math> [[fast Fourier transform]] algorithms to compute the NTT, combined with the convolution theorem, mean that the [[number-theoretic transform]] gives an efficient way to compute exact [[convolution]]s of integer sequences. While the complex DFT can perform the same task, it is susceptible to [[round-off error]] in finite-precision [[floating point]] arithmetic; the NTT has no round-off because it deals purely with fixed-size integers that can be exactly represented.

==Fast algorithms==

For the implementation of a "fast" algorithm (similar to how [[fast Fourier transform|FFT]] computes the [[discrete Fourier transform|DFT]]), it is often desirable that the transform length is also highly composite, e.g., a [[power of two]]. However, there are specialized fast Fourier transform algorithms for finite fields, such as Wang and Zhu's algorithm,<ref>[[Yao Wang]] and Xuelong Zhu, "A fast algorithm for the Fourier transform over finite fields and its VLSI implementation", IEEE Journal on Selected Areas in Communications 6(3)572–577, 1988</ref> that are efficient regardless of whether the transform length factors.

==See also==
*[[Discrete Fourier transform|Discrete Fourier transform (complex)]]
*[[Fourier transform on finite groups]]
*[[Gauss sum]]
*[[Convolution]]
*[[Least-squares spectral analysis]]
*[[Multiplication algorithm]]

==References==

<references/>

==External links==
* http://www.apfloat.org/ntt.html

{{DEFAULTSORT:Discrete Fourier Transform (General)}}
[[Category:Fourier analysis]]