Accelerometers are robust, simple to use and readily available transducers. Measuring velocity and displacement directly is not simple. In a laboratory test rig we could use one of the modern potentiometer or LVDT transducers to measure absolute displacement directly as static reference points are available. But on a moving vehicle this is not possible.
More here… OmegaArithmetic.pdf
Dear Dr Mercer:
I read your article about Omega Arithmatic and try to apply it on my experiment results, and met some problems really need your help!! I use accelerometer to collect data. Following the article, I do the Fourier transform and divided by -w^2 to get the displacement spectrum. After that when i do the inverse Fourier transform to get displacement time sequence, the result is a “complex vector”. So my question is how can I convert this “complex vector” to a displacement time sequence?? Is it correct to get the magnitude of each complex number, or just keep the real part and leave the imaginary part?? I really appreciate your kindly help!!
must do the complex calculation, use both parts. u can only use the multiplicaion and addition, basic calculation rules are:
(a+jb)*(c+jd) = a*c-b*d + j(a*d+b*c); (a+jb) + (c+jd) = (a+c) + j(b+d);
if use c, u can use two arrays to hold the real part and imaginary part separately, input spectrum vactors are acc_real[N], acc_imagine[N], N is the number of vectors.
in yr programm what u only need to calculate are: omega_disp_real[N], omega_disp_imagine[N].
Omega_disp_factor_real = -1.0/w^2; Omega_disp_factor_imagine = 0;
for each n
omega_disp_real[n] = acc_real[n]*Omega_disp_factor_real – acc_imagine[n]*Omega_disp_factor_imagine + acc_real[n]*Omega_disp_factor_imagine + acc_imagine[n]*Omega_disp_factor_real ;
omega_disp_imagine[n] = acc_real[n]*Omega_disp_factor_real – acc_imagine[n]*Omega_disp_factor_imagine;
omega_disp_real[n] = acc_real[n]*Omega_disp_factor_imagine + acc_imagine[n]*Omega_disp_factor_real ;
time squence magitude value use: magitude_disp[n] = sqrtf(omega_disp_real[n]^2 + omega_disp_imagine[n]^2); phase[n]=…
if use matlab, the vectors can be directly devided by -w^2, Omega_disp = acc_fft_vector./(-w^2).
matlab takes care of the complex calculation automatically.
time squence magitude value use magitude_disp = |Omega_disp|; phase value: …
In case of simple acceleration signal like sin(w*t), the displacement signal is displacement=IFFT(FFT(sin(w*t)/(-w^2))). It works well and good for this.
D.I (sin(w*t)) = -sin(w*t)/w^2 = IFFT(FFT(sin(w*t)/(-w^2)))
D.I. –> double integration
But however in case of composite signals with sinusoids of different frequencies and amplitudes what is the dividing factor.
D.I (sin(w*t)+sin(2*w*t)) = -(sin(w*t)/w^2)-(sin(2*w*t)/(2*w)^2)
is not equal to IFFT(FFT((sin(w*t) )+sin(2*w*t))/(-w^2)))
Sudakhar
You divide each X(f) by its own frequency f (actually divide by -[2*PI*k*df]^2)
However be careful to avoid division by zero at start – this is the dc component so set it to zero. Low frequencies are a pain as we get the so called 1/f noise – I option ignore any frequencies below 1Hz and force the to zero
Linus
The FFT and its inverse deal in complex numbers. In this sense your measured accelerations signal is the real part of a complex signal with zeroes for the imaginary part. Strictly speaking a Full FFT will give you a frequency spectrum from -half sampling rate to +half sampling rate but when we have real only data then the FFT values for the negative frequencies and the positive frequencies are related. Specifically X(-f) = conj{X(f)} that is if X(f)= a+ib then X(-f)= a-ib. In a digital representation for a real signal we have X(N-k) = conj{X(k)}.
What you get precisely when you do the FFT depends on your software package. What you need is software that has a real part only FFT which produces the usual half range result (zero to sample rate /2) and a half range Inverse FFT. This latter algorithm is often not available (it is in DATS of course!).
If you do not have the relevant inverse then
do the forward FFT of real data
divide real and imaginary parts by -(omega^2)
generate the x(N-k) complex values from conj{X(n)}
do full range inverse FFT
If done correctly you will find that the imaginary part of the inverse transform is effectively zero. You need to be very careful with the scaling. If the time between samples was dT and you have N acceleration values (N is power of two) then your frequency spacing df = 1/NdT and the value of omega for X(k) is 2*PI*k*df.
I’ve seen in training literature the statement about integration from acceleration to velocity ‘ The process involves a phase shift that supresses the higher frequency information and amplifies the lower frequency information.’ Your Omega Arithmatic article does stress the phase shift importance but does not directly confirm the statement. could you please share your comments about the relationship between phase shift and integration and the effects on the frequecy spectrum?
Integrating from acceleration to velocity, or from velocity to displacemerit, involves a 90 degree phase shift. If we integrate directly from acceleration to displacement then of course there is a 180 degree phase shift. None of these phase shifts is anything whatsoever to do with suppressing the higher frequencies. The amplitude change comes about because we are dividing by the frequency in radians/sec, usually called omega (=2 * PI * frequency). So we get the amplitude change at higher frequencies simply because the frequency is larger.
We are actually dividing by (2 * PI * i * frequency) , the ‘i’ causes the phase change and the (2 * PI * frequency) changes the amplitude. If you can bear a little maths suppose we have an acceleration spectrum A(f) with phase Theta(f). Now if R(f) is the real part and I(f) is the imaginary part then we know
A(f) = Sqrt {R(f)**2 + I(f)**2} and Theta(f) = arctan{I(f)/R(f)}
Now let us divide the read and imaginary components by (i*Omega), which is like multiplying by [1/(i*Omega)] = [-i/Omega].
Velocity Spectrum = Acceleration Spectrum * [-i/Omega]
= [R(f) + iI(f)] * [-i/Omega]
= [-i*R(f)/Omega] + [I(f)/Omega]
So the Real part is now [I(f)/Omega] and the
imaginary part is [-R(f)/Omega]
Hence the Amplitude of the Velocity Spectrum, V(f), is given by
V(f) = Sqrt{[I(f)/Omega]**2 + [-R(f)/Omega]**2}
= Sqrt {I(f)**2 + R(f)**2} / Omega
= A(f) / Omega
The velocity phase, say Phi(f), is given by
Phi(f) = arctan{[-R(f)/Omega] / [I(f)/Omega]}
= arctan{-R(f) / I(f)}
= Theta(f) – PI/2
This last relationship is not intuitive but if you do some simple geometric drawings then you can see the result.
Dear Dr. Colin Mercer
Your article Omega arithmetic is very interesting for me. I use for statistical analysis my signals (random vibration of trailer wheel from road) software Microsoft Excel (mean, standard deviation, histogram, kurtosis, skewness, FFT, …) . I measured accelerations (time domain) , I do id FFT transform (Excel can make FFT only for maximum 4096 values) and I received real and imaginary parts, then I divided each part with –omega2. I created from first 2048 values the complex values and I do it inverse FFT from 4096 parts but my Displacement has real and imaginary parts and in this point I need your help because my imaginary parts are not effectively zero. I don’t know if you know FFT in Excel.
Thank you very much for your help. I don’t know if Microsoft Excel is good software for omega arithmetic, but I will try to do it with Excel.
Excel would not be my first choice for frequency analysis! Before discussing the actual FFT in Excel, it is worth noting that doing an FFT of a random signal still has the same degree of randomness. That is there is no averaging involved.
I have looked at Excel FFT and it is a basic full range FFT.That is it computes from zero Hz to sample rate , not sample rate/2. The upper half is the complex conjugate mirror about the Nyquist frequency (sample rate/2) of the lower half. The actual frequency values you need to multiply by are for the finite half from zero to Nyquist. The frequency after to Nyquist is (Nyquist-frequency increment) and so on until you get down to (frequency increment).
For example suppose you had a sample rate of 512 samples/sec & you used a 64 point FFT then your frequency increment is 8 Hz(Sample rate / FFTsize). The frequencies then rise from 0 to 256 in steps of 8 , then as 248,240 to 8. This is because the FFF is not computing from zero to sample rate but from -Nyquist to +Nyquist. Its just the way that the FFT is implemented.
Note also that you need to divide the Excel FFT by the FFT size to correctly get to half amplitude. This can be done when you make the final time series. I do have a simple Excel Spread sheet example from displacement to velocity.
Dear Dr. Colin Mercer
Thank you very much for your answer. I sent today email to prosigusa@prosig.com with request, that my email with short example will be send directly to you.
Hi Colin,
The .pdf was very interesting, thank you for providing this resource. I have a question:
“We are actually dividing by (2 * PI * i * frequency) , the ‘i’ causes the phase change and the (2 * PI * frequency) changes the amplitude.”
Could you explain how to create this “frequency” vector? Normally, after taking a fft in Matlab, if I want to plot Velocity Amplitude vs. freq (Hz), I create a vector from 0 to #datapoints/2, then divide this vector by the length of the sample in seconds. The result is the frequency axis in (HZ).
Is a similar method used to create the frequency vector for use in Omega Arithmetic?
Thus far, my attempts to convert velocity to displacement are not working on my sine wave test signals. The phase change of the result is correct, but when I change the frequency of the velocity, the amplitude of the displacement result remains the same. Therefore I think there must be something wrong wtih my frequency vector.
Thanks for any advice,
Mark
Mark
You basically construct the frequency vector as you did for your frequency axis, together of course with the 2*pi. In principle you have to make it complex with a zero real part. Now, there is a twist, which depends on the type of FFT you do. That is, do you have a full range (0 to sample rate) or a half range FFT (0 to sample rate/2). Sample rate/2 is often called the Nyquist frequency. Also of course you want it in Real and Imaginary form.
The way you actually construct the frequency vector is similar to described in the Excel notes above.
Because of the interest l will add a Part 2 article.
Colin,
I look forward to part 2.
Matlab takes a full range FFT, as in you can see the mirror image of the data and only end up using half of it to make a regular FFT plot.
I am not having trouble constructing the positive half of the frequency vector (I don’t think) because my regular FFT plots work fine using the same method.
Per your advice I’ve tried to construct the entire frequency vector, including the negative frequencies, by mirroring my positive frequencies. But, I must have a small problem somewhere. After the omega arithmetic operation where I divide by (2*pi*freqvector*1i), most of the result’s elements are in the form xReal+Yimaginary. But, the first element is -Inf, and the last element is Inf + i*Inf.
I’ll work on it some more tomorrow…
It turns out that you can’t do an IFFT of Infinity. Changing all INF’s to 1′s got the process working.
A simple fix but I thought I’d share in case someone runs into this in the future.
Dear Dr. Colin Mercer
Thank you very much for the PDF, it was really helpful
I tried to implement whats in the PDF in Matlab to integrate a sin signal only single integration and it worked perfect match with the theoritcal integration of sin.
But when the input signal was like y = [zeros(1,16) ones(1,16)] the integrated signal wasnt as expected. the shape of the signal looks OK but the integrated values didnt look like 0000..00 12345…16 for delta T of 1 .but it was like
4.3102 4.2200 4.2699 4.2348 4.2626 ……….. 15.7626 16.7348 17.7699 18.7200 19.8102
I wonder where did this ~ 4.3102 came from? the integration should look like
0000..00 12345…16
Below is the implementation code of your idea in Matalb.
Thank you very Much !
n=32 %samples/period
delta_t = 1 %// sampling time
f=1/(delta_t*n)
periodTime = 1/f
t=0:delta_t:periodTime-delta_t; %input one period to the FFT
N=length(t)
df=1/(delta_t*N)
%this signal integration didnt work
y = [zeros(1,16) ones(1,16) ]
% below signal integration works great (commented)
%y = sin(2*pi*f*t);
% insert zero mean signal to the FFT
meanY = mean(y);
y = y – mean(y);
% construct OMEGA
for i=1:1:N/2+1
omega(i+1)=i;
end
omega(N/2+1)= 1;
for i=N/2+2:1:N
omega(i)=-(N-i+1);
end
omega(1)= 1;
omega = 2*pi*df*omega ;
trans = fft(y);
%devide each frequency by jomega
FD = trans./(j*omega);
integrated = ifft(FD) +meanY*(t+delta_t)
%theoryINT = -1*cos(2*pi*f*t)/(2*pi*f);
%plot( t,theoryINT,’-*r’,t,integrated,’-og’)
plot( t,integrated,’-O’)
Dear Dr. Colin Mercer,
Your article is very interesting and I can found some responses for my problem. In fact, I’m civil engineering and I work in a project that the principal goal it’s the road impact measurement So, to describe the impact, i measured the acceleration for each shock to obtain the road displacement. In this moment, the data processing was making in the time domain (double acceleration integration). After a baseline correction, i can obtain the signal almost correct.
I was reading about the acceleration “frequency integration” by the FFT. The problem: i work with excel. So, i need your help, I’m very interesting to obtain a copy of your excel spreadsheet to obtain the displacement by acceleration spectrum integration, because I’ve (Like Mr. Jozef Steinhubl, I get a real part and an imaginary) some problems to obtain the displacement with my spreadsheet in excel.
Thank you very much for your help.
Best regards,
Miguel Angel BENZ NAVARRETE
So, i read about the acceleration “frequency integration” by the FFT. The probleme : i work with excel. So, i nedd your help, ’cause i’m very interesting to obtain a copy of your excel spreadsheet to obtain the displacement by acceleration sprectum integration , because i’ve somes problems to obtain the displacement with my spreadsheet in excel.
Thank you very much for your help.
Dear Miguel Angel BENZ NAVARRETE
I solved the problem of FFT integration with excel. If you are interesting send me email on steinhubl.jozef@gmail.com, I can send you my excel program.
Best regards
Steinhübl
Dear Mr. Steinhubl Jozef,
I’m interesting and I will sent you an e-amil,
Je vous remercie par votre aide.
Best regards,
Miguel Angel