FFT giving undesired answer

1 view (last 30 days)
Cole Meyer
Cole Meyer on 5 Jul 2022
Edited: Star Strider on 6 Jul 2022
I'm trying to do a periodic analysis of the horizontal rows of holes on this image:
I've begun by projecting the white pixels onto the y-axis and plotting those values:
My goal is to determine the frequency of the peaks, so I went ahead with applying a fft of this data to acquire the fourier coefficients, computing the power and frequency from these (I learned from another MATLAB tutorial that power is just the magnitude of the fourier coefficients?), and plotting 1/frequency (period) against power. Finally, I believe the period corresponding to the max power should give the period that I want.
Since I've been assembling this code from examples I've seen on this forum, there are two lines I don't understand (and that I think are contributing to the issue): the lines where I calculate the maxfreq and freq. It seems that I'm trying to create a distribution of frequencies to correspond to the powers that I've calculated, but what role does maxfreq have? It seems to spread out the data in the x-direction such that I can find the maxfreq which gives any answer I want... so how do I determine what that is? Any guidance on those lines and the issue in general would be greatly appreciated!!
For reference, the code below currently gives that the vertGridSpacing is 59.94 pixels, while the true vertGridSpacing should be 123 pixels.
% Read image
bwmask = imread('bwmask.png');
% Project white pixels of bwmask onto y-axis
y_proj = sum(bwmask, 2);
x = linspace(1,length(y_proj),length(y_proj));
% Calculate fourier coefficients and eliminate first value (since it's just a sum)
tform = fft(y_proj);
tform(1) = [];
% Calculate the power corresponding to first half of fourier coefficients (since they're real)
n = length(tform);
power = abs(tform(1:floor(n/2))).^2;
% Create equally spaced frequency grid (this is the part I don't understand)
maxfreq = 1;
freq = (1:n/2)/(n/2)*maxfreq;
% Determine vertical grid spacing by taking period corresponding to maximum power
period = 1./freq;
[~, sortID] = sort(power,'descend'); % sort by x-value first
vertGridSpacing = period(sortID(1));
  6 Comments
John D'Errico
John D'Errico on 5 Jul 2022
@Cole Meyer, it looks like your real goal is to understand how to use fft to recover the period, not to find the period itself.

Sign in to comment.

Accepted Answer

Star Strider
Star Strider on 5 Jul 2022
For reference, the code below currently gives that the vertGridSpacing is 59.94 pixels, while the true vertGridSpacing should be 123 pixels.
The frequency vector of the one-sided Fourier transform should go from 0 cycles/spatial unit to the Nyquist frequency, one-half the sampling frequency (in cycles/spatial unit). So I suspect the grid spacing should be multiplied by 2 to give 119.8 instead. This is not 123 however it may be as close as the Fourier transform is able to calculate the peak.
.
  2 Comments
Star Strider
Star Strider on 5 Jul 2022
My pleasure!
I cannot run the code, since attempting to import the posted image using the online Run feature here only threw errors. One option could be to increase the frequency resolution by zero-padding the ‘yproj’ vector:
L = numel(yproj); % Assumes Vector
NFFT = 2^nextpow2(L); % For Efficiency & Increased Frequency Resolution
tform = fft(y_proj, NFFT)/L; % Slightly Changed
Fv = linspace(0, 1, NFFT/2+1)*Fn; % Frequency Vector
Iv = 1:numel(Fv); % Index Vector
promVal = 10; % Use The Appropriate Value
[pks,locs] findpeaks(abs(tform(Iv))*2, 'MinPeakProminence',promVal)
figure
plot(Fv, abs(tform(Iv)*2)
hold on
plot(Fv(locs), pks, '^r')
hold off
grid
The ‘Fn’ variable is the Nyquist frequency.
The ‘NFFT’ value can be 2 raised to the power of any positive integer multiple of the calculated nextpow2 value.
Increasing the frequency resolution is the only option I can think of that might make the frequency or period estimate of the peak value more accurate.
EDIT — (6 Jul 2022 at 5:30)
Exploring this further —
% Read image
bwmask = imread('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1056135/bwmask.png');
Warning: PNG library warning: iCCP: profile 'icc': 'RGB ': RGB color space not permitted on grayscale PNG.
% Project white pixels of bwmask onto y-axis
y_proj = sum(bwmask, 2);
Fs = 1;
Fn = Fs/2;
L = numel(y_proj); % Assumes Vector
NFFT = 2^nextpow2(L)*4; % For Efficiency & Increased Frequency Resolution
tform = fft(y_proj, NFFT)/L; % Slightly Changed
Fv = linspace(0, 1, NFFT/2+1)*Fn; % Frequency Vector
Iv = 1:numel(Fv); % Index Vector
p123 = interp1(1./Fv(2:end), abs(tform(Iv(2:end)))*2, 123);
promVal = 1E+4; % Use The Appropriate Value
[pks,locs] = findpeaks(abs(tform(Iv))*2, 'MinPeakProminence',promVal);
figure
plot(1./Fv, abs(tform(Iv))*2)
hold on
plot(1./Fv(locs), pks, '^r')
plot([1 1]*123, [0 p123], '-g')
hold off
grid
xlim([min(1./Fv) 325])
xlabel('Period (pixel)')
ylabel('Amplitude')
text(1./Fv(locs), pks, compose('\\leftarrow Period = %.1f, Amplitude = %.1f',1./Fv(locs)',pks), 'Horiz','left', 'Vert','middle')
Using:
NFFT = 2^nextpow2(L);
the period of the maximum peak (period = 128) is very close to the desired value (green vertical lline), with a difference of 5 pixels (4.065%).
Increasing that frequency resolution by 4 using:
NFFT = 2^nextpow2(L)*4;
results in the maximum peak (period = 124.1) with a difference of about 1 pixel (0.89%). It might be possible to improve that accuracy by increasing the frequency resolution even more.
The Fourier transform approach works, with longer vectors (or equivalently increasing frequency resolution by zero-padding the fft) increasing the accuracy of the estimate.
.

Sign in to comment.

More Answers (1)

Paul
Paul on 6 Jul 2022
Hi Cole,
As you've already accepted an answer, perhaps you already have what you need. But I thought it might be helpful to add this because the results that I get aren't quite the same as what you get.
Read in the image, I assume the warning doesn't affect anything.
bwmask = imread('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1056135/bwmask.png');
Warning: PNG library warning: iCCP: profile 'icc': 'RGB ': RGB color space not permitted on grayscale PNG.
% Project white pixels of bwmask onto y-axis
y_proj = sum(bwmask, 2);
figure
plot(0:(numel(y_proj)-1),y_proj)
Take the DFT
Y = fft(y_proj);
The corresponding spatial frquency vector is, assuming a sampling period of 1 pixel/sample
f = (0:(numel(Y)-1))/numel(Y); % samples/pixel
Plot the magnitude of the DFT vs. 1/f and add a line at the expected period of 123 pixels
figure
stem(1./f,abs(Y))
xline(123,'r')
xlim([100 140])
The plot shows a peak at 120 pixels (not 119.8?) and suggests an actual peak maybe between 105-140 pixels, but the spatial sampling and number of points is enough to narrow it down more than that.
As suggested by @Star Strider, zero padding can be used to fill in some of the frequencies.
Y = fft(y_proj,8192);
f = (0:(numel(Y)-1))/numel(Y); % samples/pixel
figure
stem(1./f,abs(Y))
xline(123,'r')
xlim([115 130])
Now we see the DFT "closing in" on 123 pixels, but at the end of the day you'll still have to decide how to pick the best single period for what you're trying to do. After all, as you, and @John D'Errico, have already shown, the peaks aren't equally spaced, if only by a few pixels, so there really isn't a single "period."

Community Treasure Hunt

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

Start Hunting!