How to get spatial frequency from FFT?

Hi,

I have got the first graph based on the following code. How can I get the second graph after performing FFT?

I1=0.7;

I2=0.5;

I3=0.3;

L1=200;

L2=170;

n1=1;

n2=1.444;

lam=(1.52:0.0001:1.56);

Q12=(4*pi*n1*L1)./lam;

Q23=(4*pi*n2*L2)./lam;

Q13=Q12+Q23;

I=I1+I2+I3+2*sqrt(I1*I2).*cos(Q12)+2*sqrt(I2*I3).*cos(Q23)+2*sqrt(I1*I3).*cos(Q13);

plot(lam*1000,I)

Answer 1

[Star Strider]

The Fourier transform neither knows nor cares whether the units of the independent variable are time, space, or anything else. It will do whatever you ask it to do (within limits, of course)

Try this:

I1=0.7;
I2=0.5;
I3=0.3;
L1=200;
L2=170;
n1=1;
n2=1.444;
lam=(1.52:0.0001:1.56);
Q12=(4*pi*n1*L1)./lam;
Q23=(4*pi*n2*L2)./lam;
Q13=Q12+Q23;
I=I1+I2+I3+2*sqrt(I1*I2).*cos(Q12)+2*sqrt(I2*I3).*cos(Q23)+2*sqrt(I1*I3).*cos(Q13);

figure
plot(lam*1000,I)

L = numel(lam);
Ts = mean(diff(lam));
Fs = 1/Ts;
Fn = Fs/2;
FTI = fft(I)/L;
Fv = linspace(0, 1, fix(L/2)+1)*Fn * 1E-3;
Iv = 1:numel (Fv);
[pks,locs] = findpeaks(abs(FTI(Iv)));

figure
plot(Fv, abs(FTI(Iv)))
xlim([0  0.5])
xlabel('Spatial Frequency (nm^{-1})')
ylabel('Amplitude')
text(Fv(locs), abs(FTI(locs)), sprintfc('Peak %d',(1:numel(locs))), 'HorizontalAlignment','center', 'VerticalAlignment','bottom')

producing:

.

[Sohel Rana]

Thank you so much for your help. It really worked for me. Could you please explain the following line of the code?

Fv = linspace(0, 1, fix(L/2)+1)*Fn * 1E-3;

Why multiplying by 1E-3 and Fn?

Thanks,

Sohel

[Star Strider]

My pleasure!

Sure!

The ‘Fv’ assignment creates the frequency vector for the plot. The frequency for the positive half of the symmetrical Fourier transform is defined as extending from 0 to the Nyquist frequency (‘Fn’, half the sampling frequency), and since the original linspace vector goes from 0 to 1, multiplying it by ‘Fn’ extends it from 0 to the Nyquist frequency.

The 1E-3 factor creates the correct range for the frequency axis, so it corresponds with the plot you posted. (I have no idea what the original independent variable units are in your code.)

[Star Strider]

If my Answer helped you solve your problem, please Accept it!

.

[Sohel Rana]

Hi Star,

Thank you for your feedback. Based on the code that you have generated, how can I seperate the two signals. Suppose, the first signal has a center frequency corresponding to the first peak and second signal has the center frequency corresponding to the second peak. Please note that the x-axis for the generated two signals will be the wavelength (in my code it is lam). Could you please help me with this?

[Star Strider]

As always, my pleasure!

If you have R2018a or later, use the bandpass function to create an appropriate bandpass filter for your signal. (It will be necessary to experiment to get the result you want.)

The frequencies identified by ‘Fv(locs)’ will be the centre frequencies for the peaks that findpeaks has identified.

The frequency units are irrelevant, so just use the wavelength units in the design. My code also calculates the sampling frerquency, necessary for the filter design.

[Sohel Rana]

I use R2017b. How can I do it manually instead of using bandpass function?

[Star Strider]

Since you have the Signal Processing Toolbox, make the appropriate changes to this code for each filter you want to create and use:

Fs = 44000;                                                         % Sampling Frequency
Fn = Fs/2;                                                          % Nyquist Frequency
Wp = [898 899]/Fn;                                                  % Stopband Frequency (Normalised)
Ws = [0.98 1.02].*Wp;                                               % Passband Frequency (Normalised)
Rp =  1;                                                            % Passband Ripple
Rs = 90;                                                            % Passband Ripple (Attenuation)
[n,Wp] = ellipord(Wp,Ws,Rp,Rs);                                     % Elliptic Order Calculation
[z,p,k] = ellip(n,Rp,Rs,Wp);                                        % Elliptic Filter Design: Zero-Pole-Gain 
[sos,g] = zp2sos(z,p,k);                                            % Second-Order Section For Stability

figure
freqz(sos, 2^18, Fs)                                                % Filter Bode Plot
set(subplot(2,1,1), 'XLim',Ws*Fn)                                   % Optional
set(subplot(2,1,2), 'XLim',Ws*Fn)                                   % Optional

signal_filt = filtfilt(sos, g, signal);                             % Filter Signal

You will need to create a different filter (using this code for each one separately) for each frequency band, then use it with the original signal to get the desired result.

All that is necessary is to define the passband frequencies, inside the square brackets [], in the ‘Wp’ assignment, and specify the appropriate sampling frequency ‘Fs’. This code will design an appropriate bandpass filter by default.

All the frequencies are in the same units. The code normalises them appropriately so the functions will work. If you modify it, remember that for a bandpass filter, the stopband frequencies ‘Ws’ need to include the limits of the passband frequencies, ‘Wp’. This code calculates them appropriately automatically.

Use the filtfilt function to do the actual filtering, as I do here. Supply your vector as the ‘signal’ argument.

[Sohel Rana]

Thank you very much for your detailed explanations. Actually, I wanted to the following third (center frequency corresponds to the first peak ) and fourth graphs (center frequency corresponds to the second peak )) from the first and second graphs. The x-axis would be wavelength (in my code it is lam). The first two graphs we have already got from the code and I need to seperate the last two graphs using the first two graphs since they have different center frequency.

[Star Strider]

As always, my pleasure!

You need to choose the passband frequencies from the centre frequency and the desirecd passband limits.

I am not certain what you are doing, so I must leave the details to you.

Remember the 1E-3 scaling factor in the ‘Fv’ vector in choosing the passbands with respect to your computed signals. It may be necessary to correct for that in the filter design.

[Chueng]

hello，can you explain how wavelength is converted to spatial frequency during the FFT processing? Do you have any relevant formulas?

[Star Strider]

The Fourier transform converts time, distance, or other variables to frequency units of cycles-per-original uint. So for time domain signals with the sampling frequency in seconds, the resulting frequency uints are cycles-per-second, or Hertz (Hz). Here, with the original units being nanometres, the resulting frequency uints are in cycles-per-nanometre, or more simply, or .

EDIT — (3 Jul 2024 at 14:25)

In terms of your other question (How to use fft to analyse the refelction specturm? that I just now saw), simply replace ‘cycles’ with ‘wavenumber’ since that is how you choose to express it, instead labelling it .