How to interpret a 1D-DFT of an image with a sinus grating/gradient compared to its 2D-DFT outcome?

Question

For getting a better and more intuitive understanding on how the 2D-DFT works I was playing around with sinus gratings in grayscale. I tried to compute the 1D-FFT first and compare it with the 2D-DFT to get a better sense of my happens.

So I would like to understand how the result of the 1D-DFT can be intrepreted and how this applies on my simple example. Then this knowledge should be transferred on how the 2D-DFT produces the shown outcome as well.

As far as my understanding goes the the 2D-DFT will do a FFT on the columns first and then on the rows.

So this is my code:

import numpy as np
import matplotlib.pyplot as plt

x = np.arange(-500, 501, 1)
X, Y = np.meshgrid(x, x)
wavelength = 100
angle = np.pi

grating = np.sin(2*np.pi*(X*np.cos(angle) + Y*np.sin(angle)) / wavelength)

plt.set_cmap("gray")
plt.subplot(131)
plt.imshow(grating)

oned_dft =  np.fft.fft(grating)
plt.subplot(132)

plt.imshow(np.log(abs(oned_dft)))


ft = np.fft.fft2(grating)
ft = np.fft.fftshift(ft)
plt.subplot(133)
plt.imshow(abs(ft))
plt.xlim([480, 520])
plt.ylim([520, 480])

plt.show()

I noticed that it's useful to start with angle = np.pi, since numpy doc on fft state that it will run the fft on the last axis which will a column wise operation and now the sine grating will go along the columns.

It gives the following result: When you make the image bigger you will see two straight lines along the edges and the values are decreasing column wise towards the edges and the middle.

When trying to use slightly bigger angle like angle = np.pi + 0.1 the following happens for plotting plt.imshow(np.log(abs(oned_dft)))

The lines along the rows increase when angle gets bigger, e.g. for angle = np.pi + 0.5:

It seems like the frequencies along the rows are producing further rotation of the points along the radial axis in the 2D-DFT.

How can I interpret these results? What is the underlying intuition here?

Answer 1

The 1D discrete Fourier transform (DFT) is defined as:

$$F(k)=\sum_{n=0}^{N-1}x[n]e^{-j2 \pi k n /N}$$

For periodic signals that have an integer number of periods in the $N$ samples of $x[n]$, a single value of $F(k)$ will be non-zero. This is because of the orthogonality property of the DFT. For simplicity, initially consider a complex sinusoid that oscillates with an angular frequency $2\pi k_0/N$ where $k_0$ is an integer (thus forcing it to have an integer number of periods).

$$x[n]=Ae^{j2\pi k_0 n/N}$$

The discrete Fourier transform will multiply this signal $x[n]$ by a bunch of complex sinusoids with different $k$ values, usually $k=[0,1,2,...,N-1]$.

$$F(k)=\sum_{n=0}^{N-1}(Ae^{j2\pi k_0 n/N})e^{-j2 \pi k n /N}$$ $$F(k)=A\sum_{n=0}^{N-1}e^{j2\pi (k_0-k)n/N}$$ From the DFT's clever construction, the orthogonality property means that $$k_0\ne k:\sum_{n=0}^{N-1}e^{j2\pi(k_0-k)n/N}=\frac{1-e^{j2\pi (k_0-k)}}{1-e^{j2\pi (k_0-k)/N}}=0$$ and $$k_0=k:\sum_{n=0}^{N-1}e^{j2\pi(k_0-k)n/N}=\sum_{n=0}^{N-1}e^{0}=N$$ Our signal only shows up in a single DFT bin where $k=k_0$! Therefore, the DFT of our complex sinusoid is $$F(k)=A\sum_{n=0}^{N-1}e^{j2\pi (k_0-k)n/N}=AN\delta(k_0-k)$$ where $\delta$ is the Kronecker delta function. If we had instead selected a real function, such as cosine, we can turn it into a sum of exponential terms by its definition.

$$x[n]=\cos(2\pi k_0n/N)=\frac{e^{j2\pi k_0n/N}+e^{-j2\pi k_0n/N}}{2}$$

Note that we now have two terms scaled by $\frac{1}{2}$, one with a $+k_0$ and one with a $-k_0$. By the same principle as above, we can then state

$$F(k)=A\sum_{n=0}^{N-1}\left(\frac{e^{j2\pi k_0n/N}+e^{-j2\pi k_0n/N}}{2}\right)e^{-j2\pi kn/N}=\frac{AN\delta(k_0-k)+AN\delta(k_0+k)}{2}$$

which places two peaks in the spectrum that are (conjugate) symmetric in the positive and negative frequencies.

To recap: Our DFT would only find a peak when the frequency $2\pi k/N$ matched our signal's frequency $2\pi k_0/N$.

For the 2D-DFT, we have the definition $$F(k,l)=\sum_{m=0}^{M-1}\sum_{n=0}^{N-1}x[m,n]e^{-j2\pi(km/M+ln/N)}$$

As per your problem, lets select a complex sinusoidal surface with a oscillation directed along some angle $\theta$ relative to the axis indexed by m. I'll also use $\exp$ instead of $e$ to make the exponent more clear.

$$x[m,n]=A\exp\left(j2\pi\left[k_0\cos(\theta)\frac{m}{M} + k_0\sin(\theta)\frac{n}{N}\right]\right)$$

Plugging this signal into the 2D DFT:

$$F(k,l)=A\sum_{m=0}^{M-1}\sum_{n=0}^{N-1}\exp\left(j2\pi\left[\left(k_0\cos(\theta)-k\right)\frac{m}{M} + \left(k_0\sin(\theta)-l\right)\frac{n}{N}\right]\right)$$

We can also break up the exponential terms:

$$F(k,l)=A\sum_{m=0}^{M-1}\sum_{n=0}^{N-1}\exp\left(j2\pi\left(k_0\cos(\theta)-k\right)\frac{m}{M}\right)\exp\left( j2\pi\left(k_0\sin(\theta)-l\right)\frac{n}{N}\right)$$

Assuming we carefully selected our $k_0\cos(\theta)$ and $k_0\sin(\theta)$ to be integers, we see the same as above. Applying orthogonality on each exponent turns the problem into two 1D-DFT products. The result is delta functions when $k_0\cos(\theta)-k=0$ and $k_0\sin(\theta)-l=0$. By applying orthogonality, we are left with a single peak at $F(k_0\cos(\theta),k_0\sin(\theta))$ and everywhere else is zero.

$$F(k,l)=A\sum_{m=0}^{M-1}\sum_{n=0}^{N-1}\exp\left(j2\pi\left(k_0\cos(\theta)-k\right)\frac{m}{M}\right)\exp\left( j2\pi\left(k_0\sin(\theta)-l\right)\frac{n}{N}\right)=AMN\delta(k_0\cos(\theta)-k)\delta(k_0\sin(\theta)-l)$$

Therefore, if you have shifted the DC component into the center of your 2D DFT result via fftshift, you should see a peak $k_0$ distance away from the center and at an angle $\theta$ from the x-axis. Keep in mind, imshow flips the Y-axis such that the zeroth/first element of the Y-coordinate is in the top left corner.

If your signal was real, as in the cosine example above, your result would be two (half scaled) peaks $k_0$ distance away from the center and at angles $\theta$ and $\theta+\pi$ from the x-axis. This is what you are seeing in your rightmost pictures in your question (though, zoomed in quite close).

What if the frequency you selected for your signal didn't have an integer number of periods in the window (i.e., $k_0$, or $k_0\cos(\theta)$ and $k_0\sin(\theta)$ were not integers)? You get spectral leakage where every frequency bin of the spectrum is the result of a product of a sinc function and your signal. See the plots in this post as an example. To limit the spectral leakage (if you know that there can be any significant signals with a non-integer number of periods in any dimension of your signal, you can apply a window function to your data, which diminishes the effect of spectral leakage at the tradeoff of less accurate amplitude and phase estimates of your signal of interest.