MATLAB Answers

plotting frequency distribution of fft2 transformed image...

31 views (last 30 days)
Philip on 17 Feb 2011
Answered: Edgar Guevara on 5 May 2020
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,


Show 1 older comment
Philip on 17 Feb 2011
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
David Young on 17 Feb 2011
I think you may just want to plot the power spectrum - have you had a look at fftdemo?
Philip on 17 Feb 2011
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...

Sign in to comment.

Answers (1)

Edgar Guevara
Edgar Guevara on 5 May 2020
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));
originalImage = imread(filename);
% 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;
% 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);
colormap gray;
% original image
subplot(2, 3, 1);
axis square;
set(gca, 'TickDir', 'out');
title('original image');
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);
% 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);
axis square;
set(gca, 'TickDir', 'out');
title('filtered image');
% 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)');


Sign in to comment.

Translated by