# Getting wrong frequency from fft compared to curve fit

11 views (last 30 days)
lior fridman on 31 Dec 2022
Commented: Star Strider on 31 Dec 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!

Star Strider on 31 Dec 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)
Frequencies = 1×12
-299.9891 -199.9677 -100.0215 -50.0107 -4.9709 -0.9791 0.9791 4.9709 50.0107 100.0215 199.9677 299.9891
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 CommentsShowHide 1 older comment
Star Strider on 31 Dec 2022
As always, my pleasure!

Paul on 31 Dec 2022
Edited: Paul on 31 Dec 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)
ans = 100
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.

### Categories

Find more on Signal Operations in Help Center and File Exchange

R2022a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!