plotting frequency distribution of fft2 transformed image...
10 views (last 30 days)
Show older comments
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
4 Comments
David Young
on 17 Feb 2011
I think you may just want to plot the power spectrum - have you had a look at fftdemo?
Answers (1)
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));
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
0 Comments
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!