0
$\begingroup$

I would like to numerically observe the rotation property of the Fourier transform. I believe that it is not possible since the rotation property is for the continuous transform and the DFT introduces an aliasing problem due to imperfect sampling of the input and also has a circular nature.

From the rotation property, for a rotation $\mathcal T$, $\mathcal T\mathcal F f = \mathcal F \mathcal T f$. Thus, if $f$ is symmetric to $\mathcal T$, i.e., $\mathcal T f = f$, then the rotation property implies that the transform of $f$ is symmetric to $\mathcal T$ as well, i.e., $\mathcal T\mathcal F f = \mathcal F f$.

Further, if $f$ is real, it's transform $\mathcal F f$ will be Hermitian, and along with the rotation property, this implies that $\mathcal F f$ will be real-valued, e.g., the Gaussian function (proof sketch below).

However, neither of these properties are apparent numerically, presumably due to the previously mentioned issues. Here, I am choosing $\mathcal T$ to be a $90^\circ$ rotation.

import numpy as np
np.random.seed(1)
f = np.random.randn(5, 5)
f = sum([np.rot90(f, k=k) for k in range(4)])
f.round(1)
Out[2]: 
array([[ 2.3, -1.5,  3. , -1.6,  2.3],
       [-1.6,  1.9, -4.1,  1.9, -1.5],
       [ 3. , -4.1, -1.3, -4.1,  3. ],
       [-1.5,  1.9, -4.1,  1.9, -1.6],
       [ 2.3, -1.6,  3. , -1.5,  2.3]])
np.fft.fftshift(np.fft.fft2(f)).round(1)
Out[3]: 
array([[-19.9-14.5j,  -2.7 -8.2j,   1.4 -4.4j,   7.7 -5.6j,  24.6 +0.j ],
       [ -2.9 -9.1j,  -1.9 +5.8j,  10.3 -7.5j,  -6.1 +0.j ,   6.9 +5.j ],
       [  1.4 -4.4j,  10.3 -7.5j,  -1.3 +0.j ,  10.3 +7.5j,   1.4 +4.4j],
       [  6.9 -5.j ,  -6.1 +0.j ,  10.3 +7.5j,  -1.9 -5.8j,  -2.9 +9.1j],
       [ 24.6 -0.j ,   7.7 +5.6j,   1.4 +4.4j,  -2.7 +8.2j, -19.9+14.5j]])

I also tried adding a large amount of 0 padding to see if either of these properties begins to become apparent on a "less" finite domain to try and mitigate the circular nature of the transform:

f = np.pad(f, (10000, 10000))
sum([np.abs(f - np.rot90(f, k=k)) for k in range(4)]).sum() # still symmetric
Out[4]: 5.329070518200751e-15
F = np.fft.fftshift(np.fft.fft2(f))
np.abs(np.rot90(F) - F).mean() # not symmetric
Out[5]: 12.827504372261783
np.abs(F.imag).mean() # imaginary part not 0
Out[6]: 4.852261405894498

Finally, I tried with the radially symmetric Gaussian function.

import numpy as np
x = np.concatenate([x[None] for x in np.meshgrid(*[np.linspace(-10, 10, 10000)] * 2)])
y = np.exp(-np.linalg.norm(x, axis=0))
np.abs(y - np.rot90(y)).sum() # symmetric
Out[2]: 2.091662961614887e-10
Y = np.fft.fftshift(np.fft.fft2(y))
np.abs(Y - np.rot90(Y)).mean() # not symmetric
Out[3]: 1.9977764482822713
np.abs(Y.imag).mean() # imaginary part not 0
Out[4]: 0.006017917416468471

Is it possible to numerically observe the desired results?

Sketch: If $h:\mathbb R^2\to\mathbb C$ is Hermitian and symmetric to rotations, then it is real-valued (i.e., imaginary part is $0$).

Since $h$ is symmetric to rotations, $h(x,y) = h(-x, -y)$, which implies that $\Im h(x,y)=\Im h(-x,-y)$, where $h=\Re h + i\Im h$.

Next, since $h$ is Hermitian, $\Im h(x,y)=-\Im h(-x,-y)$. Thus, $\Im h(-x,-y) = -\Im h(-x,-y)=0$, and so $\Im h(x,y)=0$

$\endgroup$

1 Answer 1

2
$\begingroup$

When I talk about the discrete Fourier transform (DFT) here, remember that the FFT is just an algorithm to realize the DFT -- so I'm talking about the FFT, too.

The usual way that the DFT is expressed has slightly weird indexing, and the way that is translated into matrices has really weird indexing. Undwinding this for the purposes of rotating 2D Fourier transforms is something that I find to be a combination of "well duh, that's obvious", and hours spent screaming at my computer while I get all of the 'i's dotted correctly, 't's crossed correctly, dotted 't's and crossed 'i's sorted out, etc.

The reason for this is that we usually start the indexing at the beginning of the sequence; there is no "negative" index. But when you start talking about rotating the input and output of the Fourier transform, you need to rotate it around the origin.

Remember that while the DFT is usually performed on vectors, the underlying math is performing its work on functions of time (or space) and frequency. Remember further that a length $N$ DFT does not change if you shift the time index by $N$ -- this is a direct consequence of the trigonometric identities $\cos \theta = \cos \theta + 2 \pi$ and $\sin \theta = \sin \theta + 2 \pi$.

In other words, the DFT works for $x(n)$ defined over any interval that is $N$ points long, as long as you define $x(n) = x(n + Nk)\, \forall\, k \in \mathbb I$.

Take the 1D case. When you "rotate" x(n) by 180$^\circ$, you need to reflect it around $n = 0$. This means that $x'(0) = x(0)$, $x'(1) = x(-1)$, etc. But the traditional indexing doesn't have an $x(-1)$. So instead, you need to calculate the indexing as $x'(0) = x(0)$, $x'(1) = x(N-1)$, etc.. This is where you start screaming at the screen...

This means that the vector $x = \begin{bmatrix}0 & 1 & 2 & 3\end{bmatrix}$, when rotated 180$^\circ$ for the discrete Fourier transform needs to end up as $x' = \begin{bmatrix}0 & 3 & 2 & 1\end{bmatrix}$

The case is similar, but more complicated, for 2D rotations. It just doesn't match the Numpy 'rot90' function (although you can wedge that into the scheme for $N \in \mathcal I^{\mathrm {odd}}$). I strongly suggest that you sketch out all four 90$^\circ$ rotations.

As an example, if you start with

$$X = \begin{bmatrix} 00 & 01 & 02 \\ 10 & 11 & 12 \\ 20 & 21 & 22 \end{bmatrix}$$

then a 90 degree rotation rotates most of your matrix "out of the frame":

$$X' = \begin{matrix} \begin{matrix} 02 & 12 & 22 \\ 01 & 11 & 21 \end{matrix} \\ \begin{bmatrix} 00 & 10 & 20 \\ \\ \\ \end{bmatrix} \end{matrix} $$ Which, in its screwy way, really means:

$$ X' = \begin{bmatrix} 00 & 10 & 20 \\ 02 & 12 & 22 \\ 01 & 11 & 21 \\ \end{bmatrix} $$


For the purposes of coding this, as mentioned, you can use Numpy fftshift if $N$ is odd, but you need to shift, then rotate, then shift back -- and it's not going to work for you for even $N$ -- for that, you need to make your own rotation function, or find one in the bowels of Numpy.

To demonstrate this:

>>> X = np.arange(25).reshape(5, 5).astype(float); X -= X.mean()
>>> X
array([[-12., -11., -10.,  -9.,  -8.],
       [ -7.,  -6.,  -5.,  -4.,  -3.],
       [ -2.,  -1.,   0.,   1.,   2.],
       [  3.,   4.,   5.,   6.,   7.],
       [  8.,   9.,  10.,  11.,  12.]])
>>> Xs = np.fft.fft2(X)
>>> Xs.round(1)
array([[  0.  +0.j , -12.5+17.2j, -12.5 +4.1j, -12.5 -4.1j, -12.5-17.2j],
       [-62.5+86.j ,   0.  -0.j ,   0.  +0.j ,   0.  +0.j ,   0.  -0.j ],
       [-62.5+20.3j,  -0.  +0.j ,  -0.  +0.j ,  -0.  +0.j ,  -0.  -0.j ],
       [-62.5-20.3j,  -0.  +0.j ,  -0.  -0.j ,  -0.  -0.j ,  -0.  -0.j ],
       [-62.5-86.j ,   0.  +0.j ,   0.  -0.j ,   0.  -0.j ,   0.  +0.j ]])
>>> np.fft.fft2(np.fft.ifftshift(np.rot90(np.fft.fftshift(X), 1))).round(1)
array([[  0.  +0.j , -62.5+86.j , -62.5+20.3j, -62.5-20.3j, -62.5-86.j ],
       [-12.5-17.2j,   0.  +0.j ,   0.  +0.j ,   0.  +0.j ,   0.  -0.j ],
       [-12.5 -4.1j,  -0.  +0.j ,  -0.  -0.j ,  -0.  +0.j ,  -0.  -0.j ],
       [-12.5 +4.1j,  -0.  +0.j ,  -0.  -0.j ,  -0.  +0.j ,  -0.  -0.j ],
       [-12.5+17.2j,   0.  +0.j ,   0.  +0.j ,   0.  +0.j ,   0.  -0.j ]])
>>> np.fft.fft2(np.fft.ifftshift(np.rot90(np.fft.fftshift(X), 2))).round(1)
array([[  0.  +0.j , -12.5-17.2j, -12.5 -4.1j, -12.5 +4.1j, -12.5+17.2j],
       [-62.5-86.j ,   0.  +0.j ,   0.  -0.j ,   0.  -0.j ,   0.  +0.j ],
       [-62.5-20.3j,  -0.  -0.j ,  -0.  -0.j ,  -0.  -0.j ,  -0.  +0.j ],
       [-62.5+20.3j,  -0.  -0.j ,  -0.  +0.j ,  -0.  +0.j ,  -0.  +0.j ],
       [-62.5+86.j ,   0.  -0.j ,   0.  +0.j ,   0.  +0.j ,   0.  -0.j ]])
$\endgroup$
4
  • $\begingroup$ How about for X = np.concatenate([x[None] for x in np.meshgrid(*[np.linspace(-2, 2, 5)] * 2)]); X = np.exp(-np.linalg.norm(X, axis=0))? Since this is a rotation symmetric function and the sample has odd dimensions, I would expect that Xs = np.fft.fft2(X); np.fft.fftshift(Xs) would also be symmetric to rotations. Further, I know that the continuous transform of this Gaussian has an imaginary part equal to 0. Any ideas why we cannot observe either here? $\endgroup$ Dec 29, 2022 at 19:53
  • $\begingroup$ I'm not making out what that code means -- sorry. Is it symmetric around $X_{1,1}$? $\endgroup$
    – TimWescott
    Dec 29, 2022 at 20:06
  • $\begingroup$ np.concatenate([x[None] for x in np.meshgrid(*[np.linspace(-2, 2, 5)] * 2)]) just creates Cartesian coordinates $\{(x,y):x,y\in\{-2,-1,0,1,2\}\}$. Then, np.exp(-np.linalg.norm(X, axis=0)) gives the Gaussian function on these coordinates, which is symmetric about 0,0 (or the $(2,2)$ entry using matrix indexing) $\endgroup$ Dec 29, 2022 at 20:20
  • 1
    $\begingroup$ OK -- you're not symmetric around the right spot. Your resulting vector of indexes should be [0, 1, 2, -2, -1]. The upper left corner of the matrix is zero-frequency (or zero position). That's what you have to index off of, or rotate around. $\endgroup$
    – TimWescott
    Dec 29, 2022 at 21:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.