Getting wrong frequency from fft compared to curve fit

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!

Answer 1

[Star Strider]

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.

.

[lior fridman]

Thank you for the detailed answer!

[Star Strider]

As always, my pleasure!

Answer 2

[Paul]

Hi @lior fridman

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.