I am trying to get the main frequency of the Fourier Transform of the image below, which is the intereference pattern detected by a CMOS sensor.

enter image description here

When I apply fft2, I got the DC peak at 0 frequency, that I know it can be suppressed if I do FT = fft2(image-mean2(image)). However, I only obtain a 0 in the (1,1) position, but the "neighbor" frequencies in the first row and column have still a higher magnitude that the peak I want to detect, which is in the (6,11) position. enter image description here

Is there a formal way to filter the image, meaning without imposing FT(1,:) = 0, FT(:,1) = 0, in order to be able to extract directly that peak?

  • $\begingroup$ Hi @Gianluca, can you post the exact image you are using? $\endgroup$
    – Bob
    Mar 12 at 10:56
  • $\begingroup$ Hey! I've just edited the post with the original image at 1080x1440, I hope it's fine now. $\endgroup$
    – Gianluca
    Mar 12 at 11:06

When processing noisy data, it is usually beneficial to filter high frequencies, this can be done with a small kernel as shown in answer by Dan Boschen.

Also if your data has an envelope it will give you some high intensity low frequencies, the low frequencies require a long kernel to be filtered, and it is better to do in the frequency domain.

The topic of filter design has been extensively explored in the past but for 1D filtering (time domain). I use the 2D information to construct a direction independent filter. Also, since the idea of this filtering is just to emphasize certain frequencies and we don't need to implement the transfer function in any recurrence equation, I just multiply the magnitude of the filter leaving the phase unchanged (I know your problem).

img = img - mean(img(:));
IMG = fft2(img);
[M, N] = size(IMG);
figure, imagesc(abs(IMG(1:64, 1:64)));

The frequencies you are interesteed in are problably those in the second spot close to 10 in the vertical axis. But we have a some frequencies that are detected higher represented at the top left corner of the image. Those need to be attenuated so that you can capture the frequencies of interest.

enter image description here

I will create a transfer function for a wide bandpass filter by cascading a a high pass filter with cutoff in 10/v, and a lowpass filter at 200/v The low frequencies are due to the envelope

f1 = 10/(N+M);;
f2 = 200/(N+M);
H = @(f) abs(f1 ./ (f1 + 1j*f) .* (1j.*f) ./ (1j*f + f2));
% see the filter response
x = linspace(0, 0.5, 1000);
figure, loglog(x, H(x));

filter response

Assuming each pixel corresponds to a square in the space (a rectanble with width = height) we can compute a frequency measured in cycles per pixel (where pixel corresponds the physical distance that separates two pixels)

% Spatial frequencies measured in cycles per pixel in each direction
phi_y = angle(exp(2j*pi*(0:M-1)/M))/(2*pi);
phi_x = angle(exp(2j*pi*(0:N-1)/N))/(2*pi);

% Construct a matrix of freq
F = sqrt(phi_y'.^2 + phi_x.^2);

If you have the actual length of each pixel you can simply divide phi_x and phi_y by the corresponding lengths.

As I said we have two options, one is to apply the analytic filter, but that changes the phase, or simply apply a real filter, by multiplying the magnitude.

Let's double check if this filtering does not modify too much the images by plotting them.

imagesc(real(ifft2(IMG .* H(F))));
colormap('gray'), colorbar();
title('Analytic filter');

imagesc(imag(ifft2(IMG .* H(F))));
colormap('gray'), colorbar();

imagesc(real(ifft2(IMG .* abs(H(F))))), colormap('gray'), colorbar();
title('Magnitude filter');

imagesc(imag(ifft2(IMG .* abs(H(F))))), colormap('gray'), colorbar();

real-imag parts

In my opinion the real filter is good, it does not leak the signal to the imaginary part and keeps the appearance of the image the same, the information you want is still there. If you want want you can increase f1 so that it filters the low frequencies even more, and you could for instance multipy by the square of the filter magnitude.

Now the last thing is to check if the peak we expect is the higher than the ones at the top left corner.

figure, imagesc(abs(IMG .* H(F))(1:64, 1:64)), colorbar();

enter image description here

  • $\begingroup$ nice work user12750353! $\endgroup$ Mar 12 at 22:03
  • $\begingroup$ When you say that you know my problem, are you referring to my previous question here? :) Could you please explain how to build the filter you proposed, and why the cutoff frequencies are defined in that way? $\endgroup$
    – Gianluca
    Mar 15 at 8:50

The mean is a low pass filter with a Sinc response which passes the 1,1 (DC) bin but as a filter has a null at every other frequency, thus when subtracted only the 1,1 bin is removed as noted.

What I would do is simply null the three undesired bins in each corner (12 bins total) with a mask (given the FFT is being computed anyway, this would be simplest) before searching for the bin with the maximum value.

Another option is to use a 2 dimensional DC nulling filter with a nulling lobe that exceeds 2 bins. This can be done by convolving a high pass filter Kernel with the image, using Matlab's imfilter command or Python with scipy.ndimage.convolve.

Here is a simple example of a simple 2D high pass filter kernel, $k$, that can be used in this case:

$$ k = \begin{bmatrix}-1/9 & -1/9 & -1/9\\-1/9 & 8/9 & -1/9\\-1/9 & 1/9 & 1/9\end{bmatrix}$$

Using python to filter with kernel k with image :

import scipy.ndimage as img
result = img.convolve(image, k)  

Using python to show the frequency response:

import scipy.signal as sp
import numpy as np
import scipy.fftpack at fft
import matplotlib.pyplot as plt
k = [[-1/9, -1/9, -1/9],[-1/9, 8/9, -1/9],[-1/9, -1/9, -1/9]]
resp = fft.fft2(k, (1080,1440))
plt.imshow(np.abs(resp), cmap = 'rainbow')

Below shows the frequency response of this simple high pass filter with the deep nulls in the 1,1 bin and adjacent bins. A more selective filter would require a much larger kernel size and simply masking the bins to ignore would be significantly less processing.

frequency response


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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