# FFT giving undesired answer

3 views (last 30 days)

Show older comments

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
on 5 Jul 2022

### Accepted Answer

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
on 5 Jul 2022

Edited: Star Strider
on 6 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');

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

.

### More Answers (1)

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

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

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

##### 0 Comments

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!