Python How can I integrate sinusoids in python code using FFT?

[jameslat]

Hello,
Thank you for taking time to read my post.

I have a discrete set of data points that represent an acceleration signal. I want to take the integral of this set of points twice so as to get a function which represents the position over time.
To accomplish this, I have taken the FFT of the acceleration signal . From there I can extract the complex harmonics and plot the signal with reasonable accuracy. This part I can do and it works great. However I keep running into problems whenever I try to integrate the sinusoids. Can you take a look at my below code and see where I might be going wrong?
Inputs:
accel, this is suppose to be the FFT of the signal we want to integrate twice
w, omega = 2 pi f
L, how many harmonics you want to use (example, L=10 would use from n=-10,-9...0...9,10)

#feed a set of FFT points.
#will alter the coefficients of the real and complex points so as to
#account for integrating each harmonic twice!
#ideally, you will be able to input the acceleration and the ouput
#will be the trajectory of the position
def int2(accel,w,L):
    sinusoids=[]
    L=int(L)
    period = (2*n.pi)/w
    xaxis = n.linspace(0,period,num=len(accel))
       
    position =[]
    tempMain = 0;
    for i in range(0,(len(accel)-1)):
        for j in range(-L,L):
            a_ = n.real(accel[j])
            b_ = n.imag(accel[j])
            n_ = len(accel)
            if j!=0:
                j2=j
                tempMain -= (1/(j2**2))*(a_*n.cos(j*w*xaxis[i])+b_*n.sin(j*w*xaxis[i]))/n_
            else:
                j2=1
                tempMain += ((xaxis[i]*xaxis[i])/(2))*(a_+b_)/n_
           
           
                       
        position.append(tempMain)
        tempMain = 0
    #we need to reverse the order of the list. This is bc the
    # algorithim that computes the FFT results in the order
    #being reversed. This is corrected whenever the IFFT is taken.
    position=position[::-1]   
       
    return position

 

[Mark44]

Hello,
Thank you for taking time to read my post.

I have a discrete set of data points that represent an acceleration signal. I want to take the integral of this set of points twice so as to get a function which represents the position over time.
To accomplish this, I have taken the FFT of the acceleration signal . From there I can extract the complex harmonics and plot the signal with reasonable accuracy. This part I can do and it works great. However I keep running into problems whenever I try to integrate the sinusoids.

What problems are you having? It would be helpful to know exactly what this code does as opposed to what you expect it to do.

Can you take a look at my below code and see where I might be going wrong?
Inputs:
accel, this is suppose to be the FFT of the signal we want to integrate twice
w, omega = 2 pi f
L, how many harmonics you want to use (example, L=10 would use from n=-10,-9...0...9,10)

#feed a set of FFT points.
#will alter the coefficients of the real and complex points so as to
#account for integrating each harmonic twice!
#ideally, you will be able to input the acceleration and the ouput
#will be the trajectory of the position
def int2(accel,w,L):
    sinusoids=[]
    L=int(L)
    period = (2*n.pi)/w
    xaxis = n.linspace(0,period,num=len(accel))
      
    position =[]
    tempMain = 0;
    for i in range(0,(len(accel)-1)):
        for j in range(-L,L):
            a_ = n.real(accel[j])
            b_ = n.imag(accel[j])
            n_ = len(accel)
            if j!=0:
                j2=j
                tempMain -= (1/(j2**2))*(a_*n.cos(j*w*xaxis[i])+b_*n.sin(j*w*xaxis[i]))/n_
            else:
                j2=1
                tempMain += ((xaxis[i]*xaxis[i])/(2))*(a_+b_)/n_
          
          
                      
        position.append(tempMain)
        tempMain = 0
    #we need to reverse the order of the list. This is bc the
    # algorithim that computes the FFT results in the order
    #being reversed. This is corrected whenever the IFFT is taken.
    position=position[::-1]  
      
    return position

 

[jameslat]

Hi Mark,

So the goal of this function is to accept the FFT of a discrete set of data points which represent an acceleration signal a(t) and find the position over time x(t).

The input is fft(a(t))=a(w).
I know that for ever harmonic in a(w), that the corresponding time domain signal is An(t) = (real(a(n))*cos(w*n*t) + imag(a(n))*sin(w*n*t))/(total number of points)
I know that if you integrate cos(nt) you get cos(nt)/n and if you integrate sin(nt) you get sin(nt)/n
There is a special case, n=0, where this will be a constant(ie n=0 and the sine signal cancel and the cosine goes to 1)

So I want to integrate everything twice. In this way I will go from acceleration to velocity and then from velocity to position.
Does that explain things well? Please let me know how I can further elaborate.

Thank you so much for your feedback :)

 

[Mark44]

I know that if you integrate cos(nt) you get cos(nt)/n and if you integrate sin(nt) you get sin(nt)/n

No to both. $\int \cos(nt)dt = \frac 1 n \sin(nt)$ and $\int \sin(nt) dt = -\frac 1 n \cos(nt)$, omitting the constants of integration. Could this be why things aren't working for you?
If you integrate the answers above you get $-\frac 1 {n^2} \cos(nt)$ for cos(nt) integrated twice and $-\frac 1 {n^2}\sin(nt)$ for the second one.