Hi Tommaso,
Starting with the acceleration as in the orginal code (I prefer to use row vectors)
% use row vectors
fs = 20; % sampling frequency
t = 0:1/fs:600; %t=t.'; % time vector
f1 = 0.01; % frequency of oscillation
acc = sin(2*pi*f1*t); % signal acceleration
N = length(acc);
Plot its fft
acc_fft = fft(acc); % fft of acceleration
f = linspace(0,fs/2,length(acc_fft)); %f=f.'; % freequency vector
figure
plot(f,abs(acc_fft),'-o'),grid
axis([0 0.02 0 1e2])
We see two problems here. First, the peak is showing up at f = 0.005, when it should be at f1 = 0.01. The frequency vector, f, is defined incorrectly. The final point is fs/2, when it should be fs - fs/N. Second, we see very small, but nonzero, points around the peak, which shouldn't happen if trying to analyze a pure sinusoid as is the case here. These points are a result of the acc vector being one element too long, i.e., the periodic extension of acc is not a perfect sine wave.
So we redefine acc
N = 12000;
t = (0:N-1)/fs;
acc = sin(2*pi*f1*t); % signal acceleration
acc_fft = fft(acc); % fft of acceleration
The sum of the acc samples should be zero, but it's not because of floating point computation errors.
sum(acc)
By definition, the first point in acc_fft should be the sum of acc, which theoretically should be zero.
acc_fft(1)
It's very small, so we will force it to be zero. This will be important later.
acc_fft(1) = 0;
Define the correct frquency vector that corresponds to acc_fft
f = (0:N-1)/N*fs; %f=f.'; % freequency vector
Plot the FFT.
figure
plot(f,abs(acc_fft),'-o'),grid
axis([0 0.02 0 1e2])
As expected, the peak is at f=0.01 and the other points are all very, very close to zero. At this point, we would be justified in setting all of the values of acc_fft to zero except the two at the peaks (you'll see the second peak at the far right edge of the plot if you change the axis command appropriately).
Now comes another issue. It looks like we are using discrete Fourier transforms/series to illustrate what happens with continuous Fourier transforms or continous complex Fourier series. However, the discrete Fourier transform, which is computed by fft, by convention, at least in Matlab, is computed for only positive frequencies (because the underlying discret time Fourier transform is periodic). However, the continuous transforms are not periodic and we have to be careful to distinguish between positive and negative frquencies. So, if we want to use the discrete Fourier transform as an analysis tool for the continuous case, we have to use fftshift and ifftshift to properly go back and forth between the discrete and continuous frequencies.
Define the frequency vector that corresponds to the continuous transform domain.
%omega = 2*pi*f; % omega vector
omega = (-N/2:(N/2-1))/N*2*pi*fs;
Because we are going to be dividing by omega later, replace omega = 0 with a small value
omega(omega == 0) = eps; % first term different from zero
fftshift acc_fft, divide by frequency, then ifftshift back, and thentake the ifft. Here, the element in acc_fft that is being divided by eps has been forced to zero above. If we hadn't done that, we'd be getting an incorrect, non-zero value in vel_fft(1), which would then cause another problem below when we do vel_fft(1)/eps.
Also, keep in mind that, in theory, only two elements of acc_fft are non-zero, and it's only those elements that should be divided by their corresponding values of omega. But, it's easier to divide the whole array by all of omega, which is why we have to protect against divide by zero on term that really shouldn't be divided at all.
vel_fft = ifftshift(fftshift(acc_fft)./(1i*omega)); % fft velocity dividing bt i omega
vel = ifft(vel_fft); % velocity in time domain
Compute the analytic results and overplot. Good match
vel_analytic = -1/(2*pi*f1)*cos(2*pi*f1*t);
figure
plot(t,real(vel));
hold on
plot(t,vel_analytic,'*')
legend('Velocità otenuta da FFT','Velocità analitica')
Do the same for displacement. Again, a key here is that the element in vel_fft that is being divided by eps is exactly zero. It it weren't, then disp_fft(1) we be a very large value, which would then shift the entire disp_fft up (or down) by that large value. I believe you're seeing that effect in your code.
disp_analytic = -1/(2*pi*f1)^2*sin(2*pi*f1*t); % analytical displacement
figure
plot(t,disp_analytic)
disp_fft = ifftshift(fftshift(vel_fft)./(1i*omega)); % fft of displacement dividing vel_fft by (i omega)
disp = ifft(disp_fft); % displacement in time domain
hold on
plot(t,real(disp))
figure
plot(t,real(disp)-disp_analytic)
Good match.
