Monday, 30 November 2015

fft - Picking the correct filter for accelerometer data


I am fairly new to DSP, and have done some research on possible filters for smoothing accelerometer data in python. An example of the type of data Ill be experiencing can be seen in the following image:


Accelerometer data


Essentially, I am looking for advice as to smooth this data to eventually convert it into velocity and displacement. I understand that accelerometers from mobile phones are extremely noisy.


I dont think I can use a Kalman filter at the moment because I cant get hold of the device to reference the noise produced by the data (I read that its essential to place the device flat and find the amount of noise from those readings?)


FFT has produced some interesting results. One of my attempts was to FFT the acceleration signal, then render low frequencies to have a absolute FFT value of 0. Then I used omega arithmetic and inverse FFT to gain a plot for velocity. The results were as follows:


Filtered signal


Is this a good way to go about things? I am trying to remove the overall noisy nature of the signal but obvious peaks such as at around 80 seconds need to be identified.



I have also tired using a low pass filter on the original accelerometer data, which has done a great job of smoothing it, but I'm not really sure where to go from here. Any guidance on where to go from here would be really helpful!


EDIT: A little bit of code:


for i in range(len(fz)): 
testing = (abs(Sz[i]))/Nz

if fz[i] < 0.05:
Sz[i]=0

Velfreq = []
Velfreqa = array(Velfreq)

Velfreqa = Sz/(2*pi*fz*1j)
Veltimed = ifft(Velfreqa)
real = Veltimed.real

So essentially, ive performed a FFT on my accelerometer data, giving Sz, filtered high frequencies out using a simple brick wall filter (I know its not ideal). Then ive use omega arithmetic on the FFT of the data. Also thanks very much to datageist for adding my images into my post :)



Answer



As pointed out by @JohnRobertson in Bag of Tricks for Denoising Signals While Maintaining Sharp Transitions, Total Variaton (TV) denoising is another good alternative if your signal is piece-wise constant. This may be the case for the accelerometer data, if your signal keeps varying between different plateaux.


Below is a Matlab code that performs TV denoising in such a signal. The code is based on the paper An Augmented Lagrangian Method for Total Variation Video Restoration. The parameters $\mu$ and $\rho$ have to be adjusted according to the noise level and signal characteristics.


If $y$ is the noisy signal and $x$ is the signal to be estimated, the function to be minimized is $\mu\|{x-y}\|^2+\|{Dx}\|_1$, where $D$ is the finite differences operator.


function denoise()


f = [-1*ones(1000,1);3*ones(100,1);1*ones(500,1);-2*ones(800,1);0*ones(900,1)];
plot(f);
axis([1 length(f) -4 4]);
title('Original');
g = f + .25*randn(length(f),1);
figure;
plot(g,'r');
title('Noisy');
axis([1 length(f) -4 4]);

fc = denoisetv(g,.5);
figure;
plot(fc,'g');
title('De-noised');
axis([1 length(f) -4 4]);

function f = denoisetv(g,mu)
I = length(g);
u = zeros(I,1);
y = zeros(I,1);

rho = 10;

eigD = abs(fftn([-1;1],[I 1])).^2;
for k=1:100
f = real(ifft(fft(mu*g+rho*Dt(u)-Dt(y))./(mu+rho*eigD)));
v = D(f)+(1/rho)*y;
u = max(abs(v)-1/rho,0).*sign(v);
y = y - rho*(u-D(f));
end


function y = D(x)
y = [diff(x);x(1)-x(end)];

function y = Dt(x)
y = [x(end)-x(1);-diff(x)];

Results:


enter image description here enter image description here enter image description here


No comments:

Post a Comment

readings - Appending 内 to a company name is read ない or うち?

For example, if I say マイクロソフト内のパートナーシップは強いです, is the 内 here read as うち or ない? Answer 「内」 in the form: 「Proper Noun + 内」 is always read 「ない...