plotting frequency distribution of fft2 transformed image...

Hi,

I wonder if anyone is able to help me find a MATLAB function/solution for plotting the Fourier transform frequencies on the x-axis (from low to high frequency), versus the average coefficient value. Or maybe someone can think of a better way of doing this..?

Basically, I would like to end up with a plot that visualises the strength of the frequencies at low frequencies against the strength at high frequencies.

I hope this makes sense,

Phil

[David Young]

What are you averaging over?

[Philip]

I am beginning to think I've asked a stupid question; I don't think the averaging part can work... I guess the result I am looking for is much like a histogram, where I want to see how many values are low frequency, how many are mid, and how many are high. So the x-axis running from the lowest to the highest frequencies should do this..?

[David Young]

I think you may just want to plot the power spectrum - have you had a look at fftdemo?

[Philip]

I have now - thanks for pointing me in its direction. The complete problem is that I have an image that has been enhanced using a low pass filter in the fft domain, and I now want to visualise the difference between the image before and after it was filtered. Unfortunately, I don't know the exact filter that was used, so I don't think I can use freqz2 to help me...

Answer 1

[Edgar Guevara]

I believe that the following example may be useful:

% spatial frequency (SF) filtering
clear; close all; clc;
filename = 'images/lena_std.tif';

% Read file and convert to grayscale
if ndims(imread(filename)) == 3
    originalImage = rgb2gray(imread(filename));
else
    originalImage = imread(filename);
end

% Compute the 2D FFT function
originalFFT = fftshift(fft2(originalImage));

% calculate the number of points for FFT (power of 2)
FFT_pts = 2 .^ ceil(log2(size(originalImage)));

%% Filter design
[f1,f2] = freqspace(FFT_pts,'meshgrid');
SF = sqrt(f1.^2 + f2.^2);
% Choose filter type
filterType = 'bandpass';
Hd = ones(size(f1));
cutoffFreq = 0.6;
cutoffFreq2 = 0.3;
Bandpass = Hd;
Lowpass = Hd;
Highpass = Hd;
Bandpass((SF < cutoffFreq2) | (SF > cutoffFreq)) = 0;
Lowpass(SF > cutoffFreq) = 0;
Highpass(SF < cutoffFreq) = 0;

%% Display parameters
show_log_amp = true; % flag to show log SF spectrum instead of SF spectrum
min_amp_to_show = 10 ^ -10; % small positive value to replace 0 for log SF spectrum

switch filterType
    case 'lowpass'
        filt = Lowpass;
    case 'bandpass'
        % SF-bandpass and orientation-unselective filter
        filt = Bandpass;
    case 'highpass'
        filt = Highpass;
end

% Perform the filtering operation in Fourier Domain
filteredFFT = filt .* originalFFT; % SF filtering
filteredImage = real(ifft2(ifftshift(filteredFFT))); % IFFT
filteredImage = filteredImage(1: size(originalImage, 1), 1: size(originalImage, 2));
filteredFFT = fftshift(fft2(filteredImage));
filteredFFT = abs(filteredFFT);

%%
figure(1);
colormap gray;

% original image
subplot(2, 3, 1);
imagesc(originalImage);
colorbar;
axis square;
set(gca, 'TickDir', 'out');
title('original image');
xlabel('x');
ylabel('y');

originalFFT = abs(originalFFT);
if show_log_amp
    originalFFT(originalFFT < min_amp_to_show) = min_amp_to_show; % avoid taking log 0
    originalFFT = log10(originalFFT);
    filteredFFT(filteredFFT < min_amp_to_show) = min_amp_to_show; % avoid taking log 0
    filteredFFT = log10(filteredFFT);
    
end




% spectral amplitude of the original image
subplot(2, 3, 2);
imagesc(mean(f1), mean(f2,2), originalFFT);
axis xy;
axis square;
set(gca, 'TickDir', 'out');
title('amplitude spectrum');
xlabel('fx (cyc/pix)');
ylabel('fy (cyc/pix)');

% filter in the SF domain
subplot(2, 3, 3);
imagesc(mean(f1), mean(f2,2), abs(filt));
axis xy;
axis square;
set(gca, 'TickDir', 'out');
title('filter in the SF domain');
xlabel('fx (cyc/pix)');
ylabel('fy (cyc/pix)');

% filtered image
subplot(2, 3, 4);
imagesc(filteredImage);
colorbar;
axis square;
set(gca, 'TickDir', 'out');
title('filtered image');
xlabel('x');
ylabel('y');

% spectral amplitude of the filtered image
subplot(2, 3, 5);
imagesc(mean(f1), mean(f2,2), filteredFFT);
axis xy;
axis square;
set(gca, 'TickDir', 'out');
title('amplitude spectrum');
xlabel('fx (cyc/pix)');
ylabel('fy (cyc/pix)');

% EOF