Matlab Accelerometer Data to Velocity

[MarkTsu]

Hello, this is my first post.

I've searched and seen some similar quesitons in the past, but thus far I'm still stuck, so I'll ask.

I'm trying to convert acceleration results into velocity using
"Omega Arithmetic". The Matlab code I wrote is something like this:

AccelFreqDomain=fft(AccelTimeSignal);
%n is the # of samples
%dt is the sampling rate
freq=(-n/2:n/2)/(n*dt) ;

VelocityFreqDomain=AccelFreqDomain./(2*pi*freq*1i);
VelocityTimeDomain=ifft(VelocityFreqDomain);
plot(VelocityTimeDomain);
hold on
plot(AccelerationTimeDomain);

When I plot the accel. and the velocity on the same graph, the phase shift is correct. When I change the amplitude of the input acceleration, the amplitude of the velocity output changes as expected. The problem is, if I change the input acceleration frequency, the output velocity amplitude remains the same, when I would expect it to change.

Any idea what I'm doing wrong?

Once I get the program responding logically, I'll probalby have some more questions about units.

Thanks for your time,
Mark

 

[IRstuff]

Is this for school?

TTFN

FAQ731-376
Chinese prisoner wins Nobel Peace Prize

 

[MarkTsu]

No. The application is engine vibration on a v6 engine. My company is the supplier of an accessory drive. I've taken a lot of accelerometer data "for fun" that's not related to our parts, and I'd like to do something with it.

They didn't teach anything close to this where I went to school..

 

[FeX32]

Is there a reason why you are using the freq. domain? You are fft then ifft.
You could use the time domain data directly and use a numerical integration technique to get an approx. velocity.
I remember using the built in function "cumsum" for this. I think it uses Simpsons rule for integrating.
I suggest you look into this.




Fe

 

[MarkTsu]

Thanks for the reply. The reasons I would like to use the frequency domain are outlined here:

http://blog.prosig.com/2010/12/16/m...tween-acceleration-velocity-and-displacement/

So far, I've used the FFT-> modify the freq domain signal-> IFFT back to time domain successfully to, say, delete all frequencies except those between 50 and 150hz or other types of filtering. My trouble is when I try to integrate or differentiate the data in the freq domain as suggested by Prosig.

 

[MarkTsu]

If you don't have time to read the whole blog, figure 12 and figure 13 summarize the difference in results between Omega Arithmetic Method and filter->integrate-> filter again.

It looks to me like there isn't much difference, especially if you were to manually remove the DC offset at the end.

But, I'd still say the Omega Arithmetic technique is "cleaner" and at least I'll understand what exactly I'm filtering. I plan to remove everyting below 10HZ (by setting = to 0 in the freq domain) before the IFFT. All the filtering options become overwhelming for a ME like myself, with the cutoff frequencies and rolloffs and phase shifts. I prefer to at least pretend I understand what I'm doing to the data than to just use some high pass filter that Matlab has available!

 

[electricpete]

freq=(-n/2:n/2)/(n*dt) ;

I think that's an error. I think the first item in the output of fft corresponds to 0 frequency (not a negative frequency). If you want to say there are negative frequencies, they live in the 2nd half of your data. When I am working with fft of a real function, I use the 1st half of the fft output and throw the 2nd half away since it is a mirrored verision of the first half.



=====================================
(2B)+(2B)' ?

 

[electricpete]

Here is an example demo of fft function I have stashed away in my notes:

% fft_try1.m - uses self-generated function
fs = 1000;       % sampling freq.
dT = 1/fs;       
TM=3;
t = 0:dT:TM;      
x = cos(2*pi*10*t) + cos(2*pi*sqrt(3)*10*t)

% example signal
%
data=x;
N=length(data);
Fmax=fs;
delta_f=1/TM;    % me: 1/TM
f = 0:delta_f:Fmax;
F=fft(data);
Fmag=abs(F)/(N/2)
figure
plot(f(1:round(N/10)),Fmag(1:round(N/10))),grid   %I added round
title('FFT'); xlabel('Frequency (Hz)');ylabel('Amplitude')
threshhold=0.2
% Now print out all the values above a certain threshhold
for i=1:round(size(F,2)/2)   
   if Fmag(i) >threshhold
      fprintf('index %g Frequency %g magnitude %g.\n',i,f(i),Fmag(i)) 

   end
   
end

=====================================
(2B)+(2B)' ?

 

[FeX32]

electricpete is on to something with that.
I don't think that line is correct either.

electricpete, btw, that code looks very familiar . It is almost exactly like I would have wrote it.
(I think I remember a thread where we all posted similar codes)

Mark,
They make a good point there. But, if you are worried about filtering the data prior to integrating or differentiating you can use the filter options in Matlab.
For example you can use this to use a butterworth filter
[b,a]=butter(3,50/(1000/2));
y=filter(b,a,data);

where a and b are the poly coef. determined by a function. Or you can determine them yourself. (y is the filtered data)
Eg. I found the coef. for a cutoff freq of 50Hz were:
% wc=50Hz
b=[0.0029 0.0087 0.0087 0.0029];
a=[1 -2.3741 1.9294 -0.5321];






Fe

 

[MarkTsu]

I think that's an error. I think the first item in the output of fft corresponds to 0 frequency (not a negative frequency).

Thanks Pete. You're right. I'm too used to working with the data after using "fftshift". From the Matlab website:

Y = fftshift(X) rearranges the outputs of fft, fft2, and fftn by moving the zero-frequency component to the center of the array.

I will investigate more tomorrow when I have the code in front of me, but I have a good feeling about the code working once I define the freq vector correctly for when fftshift has not been used.

 

[SomptingGuy]

The biggest hurdle when using Matlab's fft function is understanding the reordering of the outputs. You really need to create and test trivial inputs before moving on to real measurements. Personally I think "fftshif" just complicates things further.


- Steve

 

[FeX32]

I agree. Why do we want to shift the freq. when we know where the frequency components are already?
The code pete posted has it right on.

OP btw,
I just checked and I used cumtrapz (cumulative trapezoidal numerical integration) before with great success. It was for a impact experiment which I fixed an accelerometer to the structure. I then needed the velocity approximated from the data. Filtered then cumtrapz worked nicely. (I checked the results against a known velocities in the setup)




Fe

 

[MarkTsu]

So I have the frequency axis working ok after taking an FFT, without using fftshift. But, I'm still stuck. I think the I need to use the entire frequency vector (freq3) in the omega arithmetic : G=Y./(2*pi*freq3*1i);

Not just the first half. In the code below I tried to make the 2nd half of the frequency vector by mirroring the 1st half (minus the DC and Nyquist elements). Does anyone have any suggestions? I'll work on it more tomorrow.

%definitions
Fs=21;  %sampling rate
dt = 1/Fs; %
et = 3; % end of the interval
t = 0 : dt : et; % sampling range
y=2*sin(2 * 2 * pi * t);% define the signal

%plot in time domain
subplot(2, 1, 1); %
plot(t, y,'r.'); grid on % plot with grid
axis([0 et  -8 8]); % adjust scaling
xlabel('Time (s)'); % time expressed in seconds


%FFT secions
Y = fft(y); % compute Fourier transform enter
n = length(y)  %
Amp = abs(Y )/n; % absolute value and normalize

%if n is even the fft will be symmetric
%the first n/2 + 1 points will be unique and the rest are symmetrically
%redundant.  Point 1 is the DC component.  Point n/2 +1 is the nyquist
%component

%1 2 3 4 5  6  7  8
%4 3 2 1 0 -1 -2 -3

% if n is odd the nyquist component is not evaluated.  The number of unique
% points is (n+1)/2

%1 2 3 4  5  6  7  8   9
%4 3 2 2  1 -1  -2 -3  -4


%create the frequency axis
NumUniquePts=ceil((n+1)/2);
% the positive frequencies
freq=(0:NumUniquePts-1)*Fs/n;
%mirror the positive frequencies but make them negative
freq2=-1*freq(end:-1:1);
% delete the repetition of the nyquist component (for even length n)
freq2(1)=[];
%this is the complete frequency vector of positive and negative frequencies
freq3=[freq, freq2]
 freq3(n)=[];


%plot acceleration in frequency domain
subplot(2, 1, 2);
plot(freq, Amp(1 : NumUniquePts),'b.'); grid on % plot amplitude spectrum enter
xlabel('Frequency (Hz)'); % 1 Herz = number of cycles per second enter
ylabel('Acceleration Amplitude'); % amplitude as function of frequency enter



%Omega Arithmetic to convert acceleration to velocity
%uses complete frequency vector from above
G=Y./(2*pi*freq3*1i)

inversed=ifft(G);
figure
plot(t,real(inversed),'r')
% hold on
% plot(t,y,'b')
title('NOT WORKING'); % amplitude as function of time
xlabel('Time (s)'); % time expressed in seconds
ylabel('Velocity Amplitude');

 

[MarkTsu]

After brute force trial and error, I think this is working. Can someone advise on the units?

A 2hz sine wave with an amplitude of 2:

 y=2*sin(2 * 2 * pi * t)

is producing a velocity wave with amplitude 0.1625.

Here is the updated code for anyone curious to try it out.

 clear
home

%definitions

Fs=100; %sampling rate
dt = 1/Fs; %
et = 3; % end of the interval
t = 0 : dt : et; % sampling range
y=2*sin(2 * 2 * pi * t);% define the signal

%plot in time domain

subplot(2, 1, 1); %
plot(t, y,'r'); grid on % plot with grid
xlabel('Time (s)'); % time expressed in seconds

%FFT section

Y = fft(y); % compute Fourier transform enter
n = length(y) %
Amp = abs(Y )/n; % absolute value and normalize


%if n is even the fft will be symmetric
%the first n/2 + 1 points will be unique and the rest are symmetrically
%redundant. Point 1 is the DC component. Point n/2 +1 is the nyquist
%component

%1 2 3 4       5           6  7  8
%4 3 2 1      0          -1 -2 -3


% if n is odd the nyquist component is not evaluated. The number of unique
% points is (n+1)/2

%1 2 3 4 5  6  7  8  9
%5 4 3 2 1 -1 -2 -3 -4


%create the frequency axis

NumUniquePts=ceil((n+1)/2);

% the positive frequencies   

freq=(0:NumUniquePts-1)*Fs/(n-1);

%mirror the positive frequencies but make them negative then adjust some things
%to look right

freq2= -1*freq(end:-1:1);
if mod (n,2)==0
    freq2(1)=[];
end

freq3=[freq, freq2];
freq3';
size(freq3)

if mod(n,2)==1
freq3(n+1)=[];
end
if mod(n,2)==0
    freq3(n)=[];
end


%plot acceleration in frequency domain

subplot(2, 1, 2);
plot(freq, Amp(1 : NumUniquePts),'b.'); grid on % plot amplitude spectrum enter
xlabel('Frequency (Hz)'); % 1 Herz = number of cycles per second enter
ylabel('Acceleration Amplitude'); % amplitude as function of frequency enter


%Omega Arithmetic to convert acceleration to velocity
%uses complete frequency vector from above
G=Y./(2*pi*freq3*1i);
%result=[freq3',Y',G']   %use this line to look at freq3, Y, and G next to each other 


% Get rid of infinities resulting from divisions by zero to prepare for
% ifft
G(1)=[1];
if mod(n,2)==0
    G(n)=[1];
end



%go back to time domain with the velocity then plot it on same graph as
% acceleration
inversed=ifft(G);
figure
plot(t,real(inversed),'b')
hold on
plot(t,y,'r')
title('Maybe working?'); % amplitude as function of time
xlabel('Time (s)'); % time expressed in seconds
ylabel('Velocity Amplitude Blue, Accel RED');

 

[electricpete]

Can someone advise on the units?

If acceleration is in m/sec^2, and your algorithm divides by radian frequency in sec^-1, then you would get velocity is m/sec.

If on the other hand you started with acceleration in g's, then you'd end up with velocity in units equivalent to g-sec. You will have to do a little unit analysis if you want to convert that to something like inches per second or mm/sec.

=====================================
(2B)+(2B)' ?

 

[FeX32]

To add to the above:
"A 2hz sine wave with an amplitude of 2

is producing a velocity wave with amplitude 0.1625."

That's about right (+-2.1%)




Fe

 

[IRstuff]

Should it be d/dt(Amp*sin(2*pi*f*t) = Amp*2*pi*f*cos(2*pi*f*t)?

TTFN

FAQ731-376
Chinese prisoner wins Nobel Peace Prize

 

[electricpete]

Multiplying by i*w would accomplish the same thing in the frequency domain (actually a little bit more general since no phase is assumed ahead of time).

For this thread, we are integrating (vs differentiating), so we divide by i*w.

=====================================
(2B)+(2B)' ?

 

[electricpete]

irstuff, maybe I misunderstood your question to be a question.
If you were suggesting a check, that works:
acceleration = 2m/sec^2 at 2hz =>
velocity = 2m/sec^2 / (2*pi*2hz) = 0.159 m/sec

=====================================
(2B)+(2B)' ?

 

[IRstuff]

OK, I misunderstood his question due to his use of amplitude, instead of acceleration.

TTFN

FAQ731-376
Chinese prisoner wins Nobel Peace Prize