Matlab Accelerometer Data to Velocity
Matlab Accelerometer Data to Velocity
(OP)
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
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





RE: Matlab Accelerometer Data to Velocity
TTFN
FAQ731-376: Eng-Tips.com Forum Policies
Chinese prisoner wins Nobel Peace Prize
RE: Matlab Accelerometer Data to Velocity
They didn't teach anything close to this where I went to school..
RE: Matlab Accelerometer Data to Velocity
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
RE: Matlab Accelerometer Data to Velocity
htt
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.
RE: Matlab Accelerometer Data to Velocity
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!
RE: Matlab Accelerometer Data to Velocity
=====================================
(2B)+(2B)' ?
RE: Matlab Accelerometer Data to Velocity
CODE
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)' ?
RE: Matlab Accelerometer Data to Velocity
I don't think that line is correct either.
electricpete, btw, that code looks very familiar
(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
RE: Matlab Accelerometer Data to Velocity
Thanks Pete. You're right. I'm too used to working with the data after using "fftshift". From the Matlab website:
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.
RE: Matlab Accelerometer Data to Velocity
- Steve
RE: Matlab Accelerometer Data to Velocity
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
RE: Matlab Accelerometer Data to Velocity
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.
CODE
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');
RE: Matlab Accelerometer Data to Velocity
A 2hz sine wave with an amplitude of 2:
CODE
is producing a velocity wave with amplitude 0.1625.
Here is the updated code for anyone curious to try it out.
CODE
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');
RE: Matlab Accelerometer Data to Velocity
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)' ?
RE: Matlab Accelerometer Data to Velocity
"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
RE: Matlab Accelerometer Data to Velocity
TTFN
FAQ731-376: Eng-Tips.com Forum Policies
Chinese prisoner wins Nobel Peace Prize
RE: Matlab Accelerometer Data to Velocity
For this thread, we are integrating (vs differentiating), so we divide by i*w.
=====================================
(2B)+(2B)' ?
RE: Matlab Accelerometer Data to Velocity
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)' ?
RE: Matlab Accelerometer Data to Velocity
TTFN
FAQ731-376: Eng-Tips.com Forum Policies
Chinese prisoner wins Nobel Peace Prize
RE: Matlab Accelerometer Data to Velocity
To get a real number in the time domain, the FFT values for the negative frequencies must be the conjugate of their positive values.
M
--
Dr Michael F Platten
RE: Matlab Accelerometer Data to Velocity