Plotting Fourier series from fft output / how to obtain fundamental freq.

I want to plot my variable x_ using the fourier series, with the coefficients calculated by the fft frequency spectrum so I can change the number of terms in the series to get different numbers of harmonics with different precisions of "fit" to the data (not sure if I can use inverse fft for this instead). When I compare the fourier series plot with the original function they appear similar at first but y values are not close to each other on inspection, but when I use the Curve Fitting Tool in matlab, which does this automaticaly it gets a plot identical to the original function.

I think the problem is that I don't know how to obtain the "w" value in this equation, so i've just copied the one from the curve fitting tool. I know w is the fundamental frequency, is this just the frequency with the highest amplitude ?

eq from matlab doc https://uk.mathworks.com/help/curvefit/fourier.html

Thanks in advance!

Fs = 1000;            % Sampling frequency                    
T = 1/Fs;             % Sampling period       
L = 1000;             % Length of signal
t = (0:L-1)*T;        % Time vector#

x_ = 3*sin(2*pi*4*t)+sin(2*pi*6*t);
X = fft(x_);
terms=8;

X=(X/L)*2;
w=7;
f=zeros(1,L);
for x=1:L
    for k=2:terms % k = i in eq
        c1 = imag(X(k))*cos((k-1)*t(x)*w);
        c2 = real(X(k))*sin((k-1)*t(x)*w);
        f(x) =  f(x) + c1 + c2; 
    end
end
f=f+real(X(1));

Answer 1

[David Goodmanson]

Hi Morgan,

In this situation with periodic waves, the ftt is divided by L, as you have.

The fft is at heart an operation using complex variables. You can use sines and cosines, multiply amplitudes by 2, except don't multiply X(1) by 2, all that kind of stuff. Some of that is good for plotting purposes. But if you want to get back to the time domain it's far easier to use complex amplitudes and exp( i* ...).

The spacing of the frequency array is Fs/L. Since w = 2*pi*f, the w array is 2pi times the frequency array.

Because of aliasing, frequencies with index > L/2 can be expressed as negative frequencies but that was not done here.

In what follows I changed x_ to y (time domain), X to yf (freq domain) and f to yy (to compare with original y).

Fs =  1000;           % Sampling frequency                    
T = 1/Fs;             % Sampling period       
L = 1000;             % Length of signal
t = (0:L-1)*T;        % Time vector#

delf = Fs/L;
f = (0:L-1)*delf;     % fft output corresponds to zero frequency at element 1
w = 2*pi*f;

y = 3*sin(2*pi*4*t)+sin(2*pi*6*t);
yf = fft(y)/L;
yy = zeros(size(y));
for k = 1:L           % k is the array index
 yy = yy + yf(k)*exp(i*w(k)*t);
end
plot(t,y,t,yy+.2)     % add offset for visual purposes