plotting frequency distribution of fft2 transformed image...

조회 수: 35 (최근 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
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 CenterFile Exchange에서 Descriptive Statistics에 대해 자세히 알아보기

태그

Community Treasure Hunt

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

Start Hunting!

Translated by