Spatial Frequency (FFT) estimation from image intensity profiles

Hi,

I would like to estimate spatial frequency based on FFT from the attached image intensity data. The findpeaks() does the job though, but for some other reason I still need to find the frequency from FFT.

Answer 1

[Mathieu NOE]

hello and welcome back

let's try it with fft (have a doubt you will get a more accurate result though)

the results are in accordance with the other post (we obtained spatial period around 38 pixels so spatial freq = 1/38 = 0.026... (unit = 1/pixel)

here we get :

spatial_period_pixels = 36.5714

spatial_freq = 0.0273 +/- 0.0039

(frequency resolution is given by : df = Fs/nfft = 1/256 = 0.0039 pixel^-1

load('matlab.mat')
figure,
plot(s1)

% fft 
% first some high pass filtering 
[b,a] = butter(1,0.1,'high');
s1 = filtfilt(b,a,s1);

Fs = 1;
N = 256;
[X,freq]=positiveFFT(detrend(s1),N,Fs);
figure,
semilogy(freq,X)
% xlim([1e-2 1])
ylim([1e-4 1])
hold on 
[PKS,LOCS] = findpeaks(X,'MinPeakHeight',max(X)/5);
plot(freq(LOCS),PKS,'dr')
hold off

●

spatial_freq = freq(LOCS(1))  % first dominant (non DC) peak,  unit = 1/pixel

spatial_freq = 0.0273

spatial_period_pixels = 1/spatial_freq  % unit = pixel

spatial_period_pixels = 36.5714

%%% yet another fft function %%%%%%%%%%%%

function [X,freq]=positiveFFT(x,N,Fs)
%N=length(x); %get the number of points
k=0:N-1;     %create a vector from 0 to N-1
T=N/Fs;      %get the frequency interval
freq=k/T;    %create the frequency range
X=abs(fft(x))*2/N; % make up for the lack of 1/N in Matlab FFT

%only want the first half of the FFT, since it is redundant
cutOff = ceil(N/2); 

%take only the first half of the spectrum
X = X(1:cutOff);
freq = freq(1:cutOff);
end

[MechenG]

Thank you very much, again!. I am just curious to know, I was expecting that FFT will directly give number of bright spots in the images. Is that possible?

[Mathieu NOE]

hello again !

I am not sure how to understand what you mean by "FFT will directly give number of bright spots in the images"

here we have simply reduced the complexity of the problem by converting an image processing of image (2D array) to a signal processing (1D array)

you mean what we would get from a 2D FFT ? the other direction would not give us valuable info as there are no fringes in the other direction

[MechenG]

Hi,

Thanks again. My question was related to is there any way to identify number of fringes in the image linking with the FFT results.

[Mathieu NOE]

well the spatial representaion of the fringes is what you get from plot(s1) - or yes , a simplified vision of it as we said in the other direction there is no useful info

either we count the distance beteen the peaks or some specific level of that "signal" (as we did in the other post ) or we use fft to transform the data into spatial frequencies (as we did above)

so yes , finding the dominant peak in the fft spectrum is what tells you about the fringe spacial frequency (which is the inverse of the spacial frequency)

I am not sure to have succeeded in my explanation as I have the feeling to just repeat what I have done before (?)

[MechenG]

Yes I understood. Some has used gaussian fit (as shown in the attached image) to improve the spatial frequency resolution. Here f is the frequency corresponding to the highest peak and f* is the one estimate based on gaussian fit. I think implementation of this would further increase the resolution.

[Mathieu NOE]

ok so I tried this method , selecting the 3 dominant peaks obtained via fft and doing the gaussian fit

the frequency obtained is : fpeak_fit = 0.0379 (pixels^-1)

which would correspond to a period of : 1/fpeak_fit = 26.4056 pixels

load('matlab.mat')
figure,
plot(s1)

%% fft 
% first some high pass filtering 
[b,a] = butter(1,0.1,'high');
s1 = filtfilt(b,a,s1);
Fs = 1;
N = 256;
[X,freq]=positiveFFT(detrend(s1),N,Fs);
figure,
semilogy(freq,X)
ylim([1e-4 1])
hold on 
[PKS,LOCS] = findpeaks(X,'SortStr','descend','NPeaks',3);
fpeaks = freq(LOCS);
plot(fpeaks,PKS,'dr')

%% gaussian fit of the 3 dominant peaks
f = linspace(0.7*min(fpeaks),1.2*max(fpeaks),1000);
pf = my_gauss_fit(fpeaks,PKS,f);
semilogy(f,pf,'g')

[newpeak,ind] = max(pf);
fpeak_fit = f(ind)

fpeak_fit = 0.0379

plot(fpeak_fit,newpeak,'*r');
hold off

●

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function yfit = my_gauss_fit(x,y,xfit)

% Define g(x) = a*exp(-(x-mu)^2/(2*sigma^2)):
g = @(A,X) A(1)*exp(-(X-A(2)).^2/(2*A(3)^2));

%% ---Fit Data: Analitical Strategy---
    % Cut Gaussian bell data
    ymax = max(y); xnew=[]; ynew=[];
    ind = y > 0.05*ymax;
    xnew = x(ind);
    ynew = y(ind);

    % Fitting
    ylog=log(ynew); xlog=xnew; B=polyfit(xlog,ylog,2);

    % Compute Parameters
    sigma=sqrt(-1/(2*B(1))); mu=B(2)*sigma^2; a=exp(B(3)+mu^2/(2*sigma^2)); 

    % Plot fitting curve
    A=[a,mu,sigma];
    yfit = g(A,xfit);

end


%%% yet another fft function %%%%%%%%%%%%
function [X,freq]=positiveFFT(x,N,Fs)
%N=length(x); %get the number of points
k=0:N-1;     %create a vector from 0 to N-1
T=N/Fs;      %get the frequency interval
freq=k/T;    %create the frequency range
X=abs(fft(x))*2/N; % make up for the lack of 1/N in Matlab FFT
%only want the first half of the FFT, since it is redundant
cutOff = ceil(N/2); 
%take only the first half of the spectrum
X = X(1:cutOff);
freq = freq(1:cutOff);
end

[MechenG]

Thank you very much again.. so now we have a difference of around 10 pixels in spatial period estimated by two methods. I will explore this further and comment here my feedback..