# Omega method to integrate sin function

20 views (last 30 days)

Show older comments

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))

##### 0 Comments

### Accepted Answer

Paul
on 14 Apr 2023

Edited: Paul
on 15 Apr 2023

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.

##### 2 Comments

### More Answers (1)

Askic V
on 12 Apr 2023

Edited: Askic V
on 14 Apr 2023

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.

##### 2 Comments

Askic V
on 13 Apr 2023

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!