How to use fft to analyse the refelction specturm?

I wanna get spatial frequency from FFT, just like the picture. I have already got the reflection spectrum.The wavelength and corresponding intensity are saved in Excel.

aa = xlsread('C:\Users\jc\Desktop\6.27\xx.xlsx');

x = linspace(1510,1590,16001);%wavelength

y = aa(1:16001,2); %intensity

wavelength = x*1e-9;%nm

wavenumber = 1./wavelength;

wavenumberfit = linspace(wavenumber(1),wavenumber(16001),16001);

Answer 1

[Star Strider]

It would help to have ‘xx.xlsx’, however lacking it, I will synthesize something like it —

Fs = 0.1;
L = 1000;
t = linspace(0, Fs*L, Fs*L+1).'/Fs;
signal = sum(sin(2*pi*t*[0.125 0.215]/10).*[600 470],2);

figure
plot(t, signal)

xlabel('Wavelength (nm)')
ylabel('Optical Power (dB)')

[FTs1,Fv] = FFT1(signal,t);

figure
plot(Fv, abs(FTs1)*2)
xlabel('Spatial  Frequency (#/nm)')
ylabel('FFT Magnitude (a.u.)')

function [FTs1,Fv] = FFT1(s,t)
    t = t(:);
    L = numel(t);
    if size(s,2) == L
        s = s.';
    end
    Fs = 1/mean(diff(t));
    Fn = Fs/2;
    NFFT = 2^nextpow2(L);
    FTs = fft((s - mean(s)) .* hann(L).*ones(1,size(s,2)), NFFT)/sum(hann(L));
    Fv = Fs*(0:(NFFT/2))/NFFT;
%     Fv = linspace(0, 1, NFFT/2+1)*Fn;
    Iv = 1:numel(Fv);
    FTs1 = FTs(Iv,:);
end

My simulation is not exact, however it is close enough to demonstrate the approach. I wrote the ‘FFT1’ function for my own use, so that I wouldn’t have to type out that code whenever I wanted to calculate the FFT of a signal.

.

[Chueng]

Thank you very much for your answer. But your input is continuous , how do you need to modify it if it is discrete input? Also, I still don't understand how Fv is calculated? Why is spatial frequency calculated using this formula Fv = Fs*(0:(NFFT/2))/NFFT;

[Star Strider]

My pleasure!

First, it is not continuous. My synthesised data has asampling interval (‘Ts’) of 10 nm:

Ts = 1/Fs;

where ‘Fs’ is in units of samples/nm.

Second, frequency vector ‘Fv’ can be calculated (at least) two ways. It extends from 0 Hz (D-C) to the Nyquist frequency, Fs/2, the highest frequency that can be uniquely resolved in a sampled signal. Since the ‘complete’ Fourier transform as calculatted by the fft function extends the length of the original vector (it acutally goes from tthe ‘negative’ Nyquist frequency to the ‘positive’ Nyquist frequency if you use the fftshift function and create an appropriate frequency vector for it), calculating the half-length or one-sided Fourier transform requires using the ‘positive’ half of that frequency vector, going from 0 Hz (or whatever the frequency units are) to the Nyqust frequency. The calculation here uses the Fourier transform length ‘NFFT’ to calculate the frequency vector. It extends from 0 in steps of (1/NFFT) to NFFT/2, and then multiplies that by the sampling frequency, ‘Fs’ to get the frequency steps in the ‘Fv’ vector. If ‘NFFT’ is 10 and ‘Fs’ is 5, then:

NFFT = 10;
Fs = 5;
Fv = Fs*(0:(NFFT/2))/NFFT

Fv = 1x6

0 0.5000 1.0000 1.5000 2.0000 2.5000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Step_Length = 1/NFFT

Step_Length = 0.1000

Frequency_Vector_Step = Fs*Step_Length

Frequency_Vector_Step = 0.5000

Again, ‘frequency’ is in terms of cycles (or in this instance, wave number) per original time-domain units. So here, frequency would be in units of cycles/nm or wave-number/nm.

If you provide ‘xx.xlsx’ I can use it to write the appropriate code to calculate the actual result.

It would also help if you stated the MATLAB release/version you are using.

.

[Chueng]

Thank you for your reply.This is the data which extract from the SM125, the first column is the wavelength, and the second column is its intensity, and my MATLAB version is R2022b

[Star Strider]

Thank you.

This should work as writtten in R2022b.

The Code —

format longE
A1 = readmatrix('data.xlsx', 'HeaderLines',1)

A1 = 16001x4

1.0e+00 * 1.510000000000000e+03 -2.840000000000000e+01 NaN NaN 1.510010000000000e+03 -2.851000000000000e+01 NaN NaN 1.510010000000000e+03 -2.841000000000000e+01 NaN NaN 1.510020000000000e+03 -2.855000000000000e+01 NaN NaN 1.510020000000000e+03 -2.857000000000000e+01 NaN NaN 1.510030000000000e+03 -2.864000000000000e+01 NaN NaN 1.510030000000000e+03 -2.850000000000000e+01 NaN NaN 1.510040000000000e+03 -2.846000000000000e+01 NaN NaN 1.510040000000000e+03 -2.856000000000000e+01 NaN NaN 1.510050000000000e+03 -2.848000000000000e+01 NaN NaN

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

wavelength = A1(:,1);
intensity = A1(:,2);

Check_Indep_Var = [mean(diff(wavelength)) std(diff(wavelength))]

Check_Indep_Var = 1x2

5.000000000000000e-03 5.000156257324577e-03

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Fs = 1/mean(diff(wavelength))

Fs = 
   200

[intensityr,wavelengthr] = resample(intensity, wavelength, Fs);

Check_Indep_Varr = [mean(diff(wavelengthr)) std(diff(wavelengthr))]

Check_Indep_Varr = 1x2

5.000000000000000e-03 1.135994018793443e-13

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

% Fs = 0.1;
% L = 1000;
% t = linspace(0, Fs*L, Fs*L+1).'/Fs;
% signal = sum(sin(2*pi*t*[0.125 0.215]/10).*[600 470],2);

figure
plot(wavelengthr, intensityr)
xlabel('Wavelength (nm)')
% ylabel('Optical Power (dB)')
ylabel('Intensity (units)')
title('Original Signal')

[FTs1,Fv] = FFT1(intensityr,wavelengthr);

[pks,locs] = findpeaks(abs(FTs1)*2, 'MinPeakProminence', 0.25);

figure
plot(Fv, abs(FTs1)*2)
hold on
plot(Fv(locs), pks, 'rv')
hold off
xlim([0 2.5])
xlabel('Spatial  Frequency (#/nm)')
ylabel('FFT Magnitude (a.u.)')
text(Fv(locs), pks.', compose('\\leftarrow Mag = %.4f\n     Frq = %.4f',[pks Fv(locs).']), 'Horiz','left', 'Vert','middle')

figure
semilogx(Fv, abs(FTs1)*2)
hold on
plot(Fv(locs), pks, 'rv')
hold off
xlim([0.01 10])
xlabel('Spatial  Frequency (#/nm)')
ylabel('FFT Magnitude (a.u.)')
text(Fv(locs), pks.', compose('\\leftarrow Mag = %.4f\n     Frq = %.4f',[pks Fv(locs).']), 'Horiz','left', 'Vert','middle')

intensityr_filtered = lowpass(intensityr, 0.14, Fs, 'ImpulseResponse','iir');

figure
plot(wavelengthr, intensityr_filtered)
xlabel('Wavelength (nm)')
% ylabel('Optical Power (dB)')
ylabel('Intensity (units)')
title('Filtered Signal')

function [FTs1,Fv] = FFT1(s,t)
    t = t(:);
    L = numel(t);
    if size(s,2) == L
        s = s.';
    end
    Fs = 1/mean(diff(t));
    Fn = Fs/2;
    NFFT = 2^nextpow2(L);
    FTs = fft((s - mean(s)) .* hann(L).*ones(1,size(s,2)), NFFT)/sum(hann(L));
    Fv = Fs*(0:(NFFT/2))/NFFT;
%     Fv = linspace(0, 1, NFFT/2+1)*Fn;
    Iv = 1:numel(Fv);
    FTs1 = FTs(Iv,:);
end

The original ‘wavelength’ sampling intervals were too irregular to work with the fft function (and other signal processing applications, if you want to use them), so the original data needed to be resampled to provide consistent wavelength sampling intervals. I used the resampled data fof this analysis. There is essentially no signal energy above about 2.5 wavenumbers/nm, so I shortened the frequency axis in the plots to reflect that.

I added a lowpass filter to eliminate tthe high-frequency noise, to demonstrate how best to do that.

.

[Chueng]

Thank you very much for your reply. The reason of the irregular wavelength was because its value was extracting from the device. Therefore, the linspace function can be used to obtain the wavelength， I’ll try this to verify.

[Star Strider]

My pleasure!

Non-uniform sampling times (or wavelength values) are a simple fact-of-life in many applications. The problem is that the irregular sampling intervals make the signals unsuitable for any sort of signal processing (the only exceptin being the nufft function). I definitely recommend using the resample function instead of linspace (that likely will not provide the actual wavelengths corresponding to the recorded intensity values) or interp1 (that will only interpolate to the new wavelength values without correcting for spurious signals). The reason is that while resample interpolates to the ‘corrected’ wavelength values, it also uses an anti-aliasing filter to prevent spurious signals from appearing in the resampled signal. I recommend using my code (or some version of it) in your analysis for that reason.

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

My apologies for the delay in responding. Windows 11 crashed again and I had to restart this computer.

.