plotting frequency distribution of fft2 transformed image...

47 次查看（过去 30 天）
Philip
Philip 2011-2-17
回答： Edgar Guevara 2020-5-5
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 个评论
显示 2更早的评论
David Young
David Young 2011-2-17
I think you may just want to plot the power spectrum - have you had a look at fftdemo?
Philip
Philip 2011-2-17
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...

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Edgar Guevara
Edgar Guevara 2020-5-5
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

类别

在 Help Center 和 File Exchange 中查找有关 Fourier Analysis and Filtering 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by