- #1
jameslat
- 28
- 0
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)
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)
Code:
#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