59 views (last 30 days)
Show older comments
Chueng on 3 Jul 2024 at 13:25
Commented: Star Strider on 4 Jul 2024 at 14:34
Accepted Answer: Star Strider
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 Comments Show -2 older commentsHide -2 older comments
Show -2 older commentsHide -2 older comments
Sign in to comment.
Sign in to answer this question.
Accepted Answer
Star Strider on 3 Jul 2024 at 14:41
Edited: Star Strider on 3 Jul 2024 at 14:42
Open in 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 Comments Show 4 older commentsHide 4 older comments
Show 4 older commentsHide 4 older comments
Chueng on 4 Jul 2024 at 4:02
Direct link to this comment
https://in.mathworks.com/matlabcentral/answers/2134141-how-to-use-fft-to-analyse-the-refelction-specturm#comment_3202131
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 on 4 Jul 2024 at 10:52
Direct link to this comment
https://in.mathworks.com/matlabcentral/answers/2134141-how-to-use-fft-to-analyse-the-refelction-specturm#comment_3202441
Open in 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.
.
Chueng on 4 Jul 2024 at 12:02
Direct link to this comment
https://in.mathworks.com/matlabcentral/answers/2134141-how-to-use-fft-to-analyse-the-refelction-specturm#comment_3202521
- data.xlsx
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 on 4 Jul 2024 at 13:07
Direct link to this comment
https://in.mathworks.com/matlabcentral/answers/2134141-how-to-use-fft-to-analyse-the-refelction-specturm#comment_3202621
Open in 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>
<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 on 4 Jul 2024 at 13:52
Direct link to this comment
https://in.mathworks.com/matlabcentral/answers/2134141-how-to-use-fft-to-analyse-the-refelction-specturm#comment_3202676
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 on 4 Jul 2024 at 14:34
Direct link to this comment
https://in.mathworks.com/matlabcentral/answers/2134141-how-to-use-fft-to-analyse-the-refelction-specturm#comment_3202716
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.
.
Sign in to comment.
More Answers (0)
Sign in to answer this question.
See Also
Categories
Signal ProcessingDSP System ToolboxTransforms and Spectral AnalysisTransforms
Find more on Transforms in Help Center and File Exchange
Tags
- fft
- reflction spectrum
- spatial frenquency
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- Deutsch
- English
- Français
- United Kingdom(English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)
Contact your local office