# Frequency spectrum of a set of images in Matlab

I wish the compute the frequency distribution of a set of images. I basically am keen on finding out how the the frequency varies in the images.

I know I can use the 'fft2' function to compute the 2-dimensional Fourier Transform. I can do this for each image. But I am not sure how to obtain the isotropic frequency distribution of all the images. I want to find a method to compare the frequency distribution.

If I'm interpreting the question correctly, you're looking for a way to find a distribution of the frequencies in each images. It sounds like you might need some sort of 3D scatterplot.

If you have 20 images that are 256x256, how I imagine this working is,

1. Take your 2D FFT or 2D DCT (or blockwise FFT/DCT) for each image
2. Normalize all the values so you can have a good comparison between each image
3. Concatenate all the frequency domain images into a 256 x 256 x 20 matrix
4. Filter out the values that are below a certain threshold because they do not represent a significant component of your image. For example, if the magnitude of the values in the corners that represent high frequency information are very low, then removing them won't make a huge difference to your image because most of the information is elsewhere. But removing those data points create holes in your 3D plot so you can see the information that you're interested in
5. Plot the data in 3D. If you're using Matlab, you can use scatter3(). The plot will be easier to read if you apply a heatmap.

Otherwise, if you're looking for the frequency information of a single image, there's several ways to represent the data depending on your problem. If you clarify the question you have and it's application, I can edit this answer.

## Histogram

If you take the histogram of a grayscale image, you'll see a distribution of the intensity values of each pixel.

This easiest way to see this is to open up the image in a photo-editing software, like Photoshop, GIMP, Picasa and select histogram from the menu.

If you want to use Matlab, the hist() function plots a histogram similar to the one below. You can directly pass your image variable into the hist() function. But it's important to pick a good bin size, because different bin sizes can dramatically change the representation of the data. There are ways to calculate optimal bin sizes, but you can always guess and check.

## 2D FFT

Follow the Wikipedia article for the theoretical explanation, but basically the 2D FFT gives you an image of the intensity of low frequency to high frequency information. The center of the image is low frequency information (such as flat surface, sky, skin, walls) and the corners are high frequency information (edges, grainy noise, intricate patterns).

If you're using Matlab, you need to shift the 2D FFT to put the low frequency information in the center of the image. It's also a good idea to normalize the 2D FFT so you can see the details better; it's similar to contrast enhancing a photograph to see the details.

## 2D DCT

Follow the Wikipedia article for the theoretical explanation. The 2D DCT is very similar to the 2D FFT; It also gives you an image of the intensity of low frequency to high frequency information, but the low frequency information is in the top-left corner and the high frequency information is in the bottom-right corner. This makes the information easier to manipulate (at least for me).

Generally, the 2D DCT is used on blocks of an image (8x8 pixels). This gives you frequency information that is still related to its location in the spatial image.

The 2D DCT is useful for encoding images. A variant of JPEG encoding uses the 2D DCT transform.

Projection onto 1D

This is similar to what I use, it is an attempt to project the 2D spectrum onto 1D. That is, it tries to perform

$$\hat{f}_{pol}(r) = \int_{0}^{2\pi} \hat{f}(r,\theta) \, d\theta$$

where $[u,v] = [r \cos\theta, r \sin\theta]$ and $u,v$ are the spectrum coordinates. There is a section in the code where the spectrum is multiplied by $r$, this is to try to adjust for the interpolation errors close to the centre due to discretisation. The 1D projections can then be compared between images.

close all

% Image
imagesc(I); colormap gray; title('Original image'); pause;

% Fourier transform
F = ifftshift(fft2(I))./rows./cols;

% Show spectrum (log)
imagesc(log(abs(F))); title('Fourier transform (abs log)'); pause;

% Grid of FFT coordinates
[rows, cols] = size(F);
[ux, uy] = meshgrid(([1:cols]-(fix(cols/2)+1))/(cols-mod(cols,2)), ...
([1:rows]-(fix(rows/2)+1))/(rows-mod(rows,2)));

% Convert to polar coordinates
th = atan2(uy,ux);
r = sqrt(ux.^2 + uy.^2);

% Convert to polar coordinates
Fr = F .* r;
imagesc(abs(Fr)); title('Fourier transform x radius'); pause;
rcoords = linspace(0,sqrt(ux(1,1)^2 + uy(1,1)^2),rows);
thcoords = linspace(0,2*pi,cols);
[ri,thi] = meshgrid(rcoords,thcoords);
[x,y] = pol2cart(thi,ri);
Fp = interp2(ux,uy,abs(Fr),x,y);
imagesc(Fp); title('Fourier transform in polar coordinates'); pause;

% Sum columns to give 1D projection
F1D = sum(Fp);
plot(rcoords,F1D); title('Projection onto 1D'); xlim([0 0.5]);


Though I do not have experience using this particular function on images, sometimes in communications we use a spectrogram to see how the frequency distribution of a signal changes over time take a look at this wiki page http://en.wikipedia.org/wiki/Spectrogram.

Essentially on a one dimensional array of size 2048 we could take a set of smaller FFTs say size 16 of size 128 with out over lap and plot the ffts over time. Other variations include using overlap when taking the fft.

This would be difficult to use in image processing because of the fourth dimension.

However, one way you could reduce the dimensionality of something to this effect over time is by taking an FFT of a single row of the image over the entire set of the images and plot these as a three dimensional surface, then do this for each relevant row in the image.