Decomposes an image into fourier 2D fourier components and reconstructs it via 2D Fourier Series using a chosen relative bandwidth. 100% of the bandwidth recorvers the original image
Reads an image as a 2D Numpy Array, which we will handle as a (pseudo)tensor
The Discrete Fourier Transform of a sequence is given by
And the inverse discrete fourier transform is the fourier series.
Which can be written as a tensor product (using Einstein's Summation Notation)
Where
To get a fourier series using only part of the bandwidth we can just limit the frequencies k we use.
f = np.fft.fft(signal)
f = np.fft.fftshift(f)
It's worth mentioning that since we are working with complex fourier series, we are also using negative frequencies. The defined bandwidth contains the negative frequencies as well. That is the reason we are shifting
To create $E^{nk} we use a mesh. For that we generate all values of
band_slice = int(relative_bw * N)
And with that we generate our n and k arrays for the mesh
n_array = np.arange(0, N)
k_array = np.arange(-int(band_slice / 2), int(band_slice / 2))
We create the mesh
nk_mesh = np.meshgrid(n_array, k_array, indexing="ij")
And use it to create the tensor with the help of numpy. The first component of the mesh are the n-values and the second are the k-values
e_tensor = np.exp(1j * 2*np.pi / N * nk_mesh[0]*nk_mesh[1])
We don't forget to slice our used frequencies
slice_index = int((N + band_slice) / 2)
f_k = f[-slice_index : slice_index]
And now we can finally multiply the tensors
x_n = np.tensordot(f_k, e_tensor, axes=(0, 1))
x_n = x_n / N
The 2D Fourier series of a 2D Array is
As a tensor product
The problem with this is that allocating a 4D array requires Petabytes of memor, which is why the
First, we find the product
Then
and finally
The generation of the meshes and tensors is analogue to the 1D version but using tuples of meshes. The final tensor product is
x_nm = np.tensordot(f_k_tensor, exp_tensor[0], axes=(0, 1))
x_nm = np.tensordot(x_nm, exp_tensor[1], axes=(0, 1))
x_nm = x_nm / (N[0] * N[1])