FFT giving undesired answer

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));

[John D'Errico]

I always laugh whenever I see someone telling us FFT is giving the wrong answer or there is a bug in double precision arithmetic, etc.

FFT is giving the answer it is designed to return. You may be interpreting the result of what you passed in incorrectly.

[Cole Meyer]

Thank you for your response, John. You're absolute right. I understand that I am either passing the function something wrong or interpreting the result wrong – that's what I'm hoping to troubleshoot. I've changed the title of my question to reflect this.

[Paul]

Hi Cole,

Posting the file bwmask.png will allow you to get more help because then people can run your code and/or post their own solution for the same input image.

You can edit the Question and click the paper clip icon in the "Insert" tools if you want to go down this path.

BTW, if all you want to do is pull out the primary frequency (or period) there's no reason to bother with squaring the coefficients, which is actually energy, nor the one-sided stuff. Just plotting the magnitude of the coefficients vs the full frequency vector should be all that is needed for this problem.

[Cole Meyer]

Thank you for your advice, Paul – I've just attached bwmask.png and simplified my code (outside of this post) using your suggestions! As I mentioned in my reply to Star Strider, this unfortunately didn't solve my issue, but is much cleaner nonetheless.

EDIT: Perhaps it would help if I demonstrate why I'm not satisfied with the result of 119.8. As the image below shows, the peak values give spacings of 123, 123, 123, 122, 124, 122, and 125, which have a mean of 123 and a std of 0.99.

[John D'Errico]

@Paul Actually, you can just save the image shown, then attach it, as I did. Then you can use the image too, since MATLAB online can see and use attachments.

[John D'Errico]

@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.

Answer 1

[Star Strider]

‘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.

.

[Cole Meyer]

Thank you for your response! It would make much more sense that the grid spacing should be multiplied by 2 – that helps. As for the slight discrepency between 119.8 and 123, can you spot a place within my code that I might be able to increase the resolution? My goal is to fit a grid to the holes, so spacing that's off by even 3 pixels will compound with each additional row.

[Star Strider]

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.

.

Answer 2

[Paul]

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."