Omega method to integrate sin function

Hi, I'm a student who is practicing with signal processing and matlab. I'm trying to integrate a sine function dividing it by (i*2*pi*f). And I'm trying to do that two times as if my signal was an acceleration and I would like to calculate displacement. I can't understand why it works to obtain velocity but it doesn't work with second integration. This is the code. Thank very much to everyone.

clear all; close all; clc;

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(t,acc);

acc_fft = fft(acc); % fft of acceleration
f = linspace(0,fs/2,length(acc_fft)); f=f.'; % freequency vector
omega = 2*pi*f; % omega vector
omega(1) = eps; % first term different from zero
vel_fft = acc_fft./(1i*omega); % fft velocity dividing bt i omega
vel = ifft(vel_fft); % velocity in time domain
vel_analytic = -1/(2*pi*f1)*cos(2*pi*f1*t); 

plot(t,real(vel));
hold on
plot(t,vel_analytic,'*')
legend('Velocità otenuta da FFT','Velocità analitica')

disp_analytic = -1/(2*pi*f1)^2*sin(2*pi*f1*t); % analytical displacement
figure
plot(t,disp_analytic)
disp_fft = 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))

Answer 1

[Paul]

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)

ans = -2.5463e-13

By definition, the first point in acc_fft should be the sum of acc, which theoretically should be zero.

acc_fft(1)

ans = -6.4812e-15

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.

[Askic V]

Wow, excellent and detail explanation @Paul. Always so much to learn from your posts. Thank you for that!

[Tommaso Ballantini]

Thank you very much @Paul.

Answer 2

[Askic V]

You're getting the weird values for disp because of division with eps (very small number).

This code will produce what you expect:

clear all; close all; clc;
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(t,acc);

acc_fft = fft(acc); % fft of acceleration
f = linspace(0,fs/2,length(acc_fft)); f=f.'; % freequency vector
omega = 2*pi*f; % omega vector
omega(1) = 1; % first term different from zero
vel_fft = acc_fft./(1i*omega); % fft velocity dividing bt i omega
vel = ifft(vel_fft); % velocity in time domain
vel_analytic = -1/(2*pi*f1)*cos(2*pi*f1*t); 
plot(t,real(vel));
hold on
plot(t,vel_analytic,'*')
legend('Velocità otenuta da FFT','Velocità analitica')

disp_analytic = -1/(2*pi*f1)^2*sin(2*pi*f1*t); % analytical displacement
figure
plot(t,disp_analytic)
disp_fft = vel_fft./(2i*omega); % fft of displacement dividing vel_fft by (i omega)
disp = ifft(disp_fft); % displacement in time domain
hold on
plot(t,real(disp),'o')

Please note there is a division with 2*i*omega for finding disp_fft.

I'm not sure why is that, but it may have something to do with even/odd functions. There are many people here who are very familiar with FFT and IFFT and would probably know why this scaling is important.

@John D'Errico @Walter Roberson or @Paul

[Tommaso Ballantini]

Thank you very much for your reply.

P.s. my fault, there should not have been a 2.

[Askic V]

No, it is not your fault, I added disp_fft = vel_fft./(2i*omega) in order to scale it properly. Your original code had 1i*omega. Execute and see for yourself.