44 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Chueng am 3 Jul. 2024 um 13:25
Kommentiert: Star Strider vor etwa 2 Stunden
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 Kommentare -2 ältere Kommentare anzeigen-2 ältere Kommentare ausblenden
-2 ältere Kommentare anzeigen-2 ältere Kommentare ausblenden
Melden Sie sich an, um zu kommentieren.
Melden Sie sich an, um diese Frage zu beantworten.
Antworten (1)
Star Strider am 3 Jul. 2024 um 14:41
Bearbeitet: Star Strider am 3 Jul. 2024 um 14:42
In MATLAB Online öffnen
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 Kommentare 4 ältere Kommentare anzeigen4 ältere Kommentare ausblenden
4 ältere Kommentare anzeigen4 ältere Kommentare ausblenden
Chueng vor etwa 13 Stunden
Direkter Link zu diesem Kommentar
https://de.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 vor etwa 6 Stunden
Direkter Link zu diesem Kommentar
https://de.mathworks.com/matlabcentral/answers/2134141-how-to-use-fft-to-analyse-the-refelction-specturm#comment_3202441
In MATLAB Online öffnen
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 vor etwa 5 Stunden
Direkter Link zu diesem Kommentar
https://de.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 vor etwa 4 Stunden
Direkter Link zu diesem Kommentar
https://de.mathworks.com/matlabcentral/answers/2134141-how-to-use-fft-to-analyse-the-refelction-specturm#comment_3202621
In MATLAB Online öffnen
- 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 vor etwa 3 Stunden
Direkter Link zu diesem Kommentar
https://de.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 vor etwa 2 Stunden
Direkter Link zu diesem Kommentar
https://de.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.
.
Melden Sie sich an, um zu kommentieren.
Melden Sie sich an, um diese Frage zu beantworten.
Siehe auch
Kategorien
Signal ProcessingDSP System ToolboxTransforms and Spectral AnalysisTransforms
Mehr zu Transforms finden Sie in Help Center und 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!
Es ist ein Fehler aufgetreten
Da Änderungen an der Seite vorgenommen wurden, kann diese Aktion nicht abgeschlossen werden. Laden Sie die Seite neu, um sie im aktualisierten Zustand anzuzeigen.
Website auswählen
Wählen Sie eine Website aus, um übersetzte Inhalte (sofern verfügbar) sowie lokale Veranstaltungen und Angebote anzuzeigen. Auf der Grundlage Ihres Standorts empfehlen wir Ihnen die folgende Auswahl: .
Sie können auch eine Website aus der folgenden Liste auswählen:
Amerika
- América Latina (Español)
- Canada (English)
- United States (English)
Europa
- 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)
Asien-Pazifik
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)
Kontakt zu Ihrer lokalen Niederlassung