Getting wrong frequency from fft compared to curve fit
10 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
lior fridman
am 31 Dez. 2022
Kommentiert: Star Strider
am 31 Dez. 2022
I'm trying to get frequency of a signal using fft and plotting its spectrum.
However I get slightly off frequency and cannot figure out why, I've even tried to do a simple fit to my signal and got the correct frequnecy from it.
I suspect that I'm creating a wrong frequency vector which causes this deviation, but can't tell for certain.
This is the relevent code -
%%%
lambda = 1000e-9;
k = 2*pi/lambda;
delta_k = 5000;
N = 100;
k_vec = (k-(N)/2*delta_k):delta_k:(k+(N-1)/2*delta_k);
A = 0.1; % reflection coef.
delta_x = 50e-6; %m
E_mirror = 1;
E_obj = A*exp(1i*k_vec*2*delta_x);
I = abs(E_mirror+E_obj).^2;
I_mirror = abs(E_mirror)^2;
I_obj = abs(E_obj).^2;
I_diff = I-I_mirror-I_obj;
F_I_diff = fft(I_diff);
z_scale = 1/delta_k;
z_vec = z_scale/N*(-N/2:N/2-1);
%%%
The signal which frequency Im trying to find is I_diff.
Help would be appreciated, Thanks!
0 Kommentare
Akzeptierte Antwort
Star Strider
am 31 Dez. 2022
Calculating the frequencies for a two-sided Fourier transform can initially be a bit of a challenge.
Try something like this —
Fs = 1234; % Sampling Frequency
Ls = 10;
t = linspace(0, Fs*Ls-1, Fs*Ls)/Fs; % Time Vector
s = sum(sin([1; 5; 50; 100; 200; 300]*2*pi*t)); % Signal Vector
figure
plot(t, s)
grid
xlabel('Time')
ylabel('Amplitude')
title('Signal')
Fn = Fs/2; % Nyquist Frequency
L = numel(t); % Signal Length
NFFT = 2^nextpow2(L); % For Efficiency
FTs = fft(s,NFFT)/L; % Fourier Transform
FTss = fftshift(FTs); % Shift For Central Symmetry
Fv = linspace(-Fn, Fn-Fs/NFFT, NFFT); % Frequency Vector For Two-Sided 'fft' With An Even Number Of Elements
figure
plot(Fv,abs(FTss))
grid
xlabel('Frequency')
ylabel('Magnitude')
title('Fourier Transform')
[pks,locs] = findpeaks(abs(FTss), 'MinPeakProminence',0.25); % Use 'findpeaks' To Return Peak Indices To Demonstrate Peak Frequencies
Frequencies = Fv(locs)
figure
plot(Fv,abs(FTss))
grid
axis([[-1 1]*10 0 0.7])
xlabel('Frequency')
ylabel('Magnitude')
title('Fourier Transform (Detail)')
txtstr = compose('\\leftarrow %+9.4f Hz',Frequencies);
text(Frequencies(5:8), pks(5:8), txtstr(5:8), 'Vert','top', 'Horiz','left', 'Rotation',45)
Using the ‘NFFT’ value as calculated here serves two purposes. First, it improves the efficiency of the fft calculation, and second it creates a fft vector with an even number of elements.
.
2 Kommentare
Weitere Antworten (1)
Paul
am 31 Dez. 2022
Bearbeitet: Paul
am 31 Dez. 2022
Assuming you're only interested in the magnitude of the DFT, it looks like everything is correct except for not applying fftshift to the output from fft. Howver, if you're also interested in the phase, then an additional step may be required because the first point in I_diff does not correspond to k = 0.
lambda = 1000e-9;
k = 2*pi/lambda;
delta_k = 5000;
N = 100;
k_vec = (k-(N)/2*delta_k):delta_k:(k+(N-1)/2*delta_k);
A = 0.1; % reflection coef.
delta_x = 50e-6; %m
E_mirror = 1;
E_obj = A*exp(1i*k_vec*2*delta_x);
I = abs(E_mirror+E_obj).^2;
I_mirror = abs(E_mirror)^2;
I_obj = abs(E_obj).^2;
I_diff = I-I_mirror-I_obj;
figure
plot(k_vec,I_diff),ylabel('Idiff')
F_I_diff = fft(I_diff);
z_scale = 1/delta_k;
zvec is correc for an even-length, fftshifted DFT
numel(I_diff)
z_vec = z_scale/N*(-N/2:N/2-1);
F_I_diff = fftshift(F_I_diff);
plot(z_vec,abs(F_I_diff))
The peaks seem to correspond to the input signal I_diff. My by-eye estimate is that the peaks should be at roughly +- 1.5e-5.
0 Kommentare
Siehe auch
Kategorien
Mehr zu Discrete Fourier and Cosine Transforms finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!