How to use fft to analyse the refelction specturm?

Question

Chueng 2024-7-3，13:25

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2134141-how-to-use-fft-to-analyse-the-refelction-specturm

评论： Star Strider 2024-7-4，14:34

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2024-7-3，14:41

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2134141-how-to-use-fft-to-analyse-the-refelction-specturm#answer_1480681

编辑：Star Strider 2024-7-3，14:42

在 MATLAB Online 中打开

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.

.

6 个评论
显示 4更早的评论隐藏 4更早的评论

Star Strider 2024-7-4，10:52

在 MATLAB Online 中打开

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.

.

Star Strider 2024-7-4，13:07

在 MATLAB Online 中打开

data.xlsx

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>

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

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

% 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 2024-7-4，13:52

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 2024-7-4，14:34

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.

.

请先登录，再进行评论。