filtfilt changes signal amplitude badly - how to choose the right filter?

Hi all,

I have been trying to filter a the displacement traces of a moving animal to calculate kinematics such as velocity and acceleration. The problem is that I do not know anything about how the signal should appear in the time or frequency domain - the animal moves continously varrying it's speed. I tried using various filters to remove the jitter and they seem to dampen the amplitude of the velocity a lot as you can see in the image (instantaneous velocity in blue is the unfiltered velocity trace and the filtered one is in magenta). Can anyone who knows more about signal processing help me to choose the right filter which can smothen this signal without loosing much information?

I am attaching a sample code and the starting variable xPos and yPos which will bring you on the same page as me in just few seconds.

Thanks already. :)

%%%%%%%%%%%%%%%%%%%%%%%
dx   = [0; diff(xPos)]; % distance between two consecutive x-coordinates
dy   = [0; diff(yPos)]; % distance between two consecutive y-coordinates

camscale_px_per_mm = 22;
datarate_Hz=700;

tmp_dist_unfilt           = sqrt(dx.^2 + dy.^2);
tmp_dist_unfilt           = tmp_dist_unfilt./camscale_px_per_mm;  % convert to distance in mm
tmp_vel_unfilt            = tmp_dist_unfilt.*datarate_Hz;         % convert to velocity in mm/s

filterF = ones(1,4)./4; %for filtering out jitter

idx_nan     = isnan(dx);
idx_nan = find(idx_nan==1);

dx(idx_nan) = 0; % for filtering nan values need to be removed 
dy(idx_nan) = 0;

dxF        = filtfilt(filterF, 1, dx);
dyF        = filtfilt(filterF, 1, dy);

tmp_dist_fF = sqrt(dxF.^2 + dyF.^2);     % distance moved between iterations, in pixels
tmp_dist_fF = tmp_dist_fF./camscale_px_per_mm;  % convert to distance in mm
tmp_vel_fF  = tmp_dist_fF.*datarate_Hz;  % convert to velocity in mm/s

tmp_vel_fF(idx_nan) = nan; % re-insert the nan values

yyaxis left
plot(tmp_vel_unfilt);
yyaxis right
hold on;plot(tmp_vel_fF);
%%%%%%%%%%%%%%%%%%%%%%%%%

[Jan]

What are the inputs? Did you measure the positions directly? What is the wanted output? The average or locally smoothed velocity?

Of course the moving average filter you have applied damps the amplitude. To suggest a "better" filter, you have to define, what you consider as wanted signal. Afterwards you can design a filter. Without defining, what the noise is, it is impossible to remove it without destroying the wanted information.

Please explain, what you want to achieve.

If you want to get a smoother velocity, you can use gradient instead of diff, because it uses a 2nd order quotient of differences instead of a 1st order. A larger time step helps also directly:

diff_X  = X(2:end) - X(1:end-1);        % same as diff(x)
diff_X8 = (X(9:end) - X(1:end-8)) / 8

This obtains the average velocity over 8 time steps without smoothing the original signal.

[Marion Rosello]

Thanks for your comments, Jan.

The positions: I use a high speed (700Hz) camera to capture the images and do some image processing on the fly to determine the centroid of thresholded images. I think the jitter comes from this small variations in determining the right pixels of the centroid.

Output: I'm looking for a locally smoothed velocity trace.

Answer 1

[Bjorn Gustavsson]

Don't instantly start looking at your data after differentiation - differentiation is a noise-amplifying operation! First you have to take a look at the data and have a think about what your data and its intrinsic noise-characteristics:

t = (1:numel(xPos))/700;
subplot(2,2,1)
plot(t,xPos)
subplot(2,2,2)
plot(t,yPos)
subplot(2,2,3)
plot(xPos,yPos)
subplot(2,2,4)
comet(xPos,yPos)

Then you'll have to make some kind of decision on whether your animal actually moves in 1/700th of a second or if positional variations at that time-scale is just noise (Just look at the variation in kinetic energy or try to calculate the work done by your animal - integrate scalar product of velocity and acceleration divided by mass, is that even plausible considering the metabolism of you subject? Are the accelerations consistent with the maximum/steady force your animal can generate).

I'd guess you should filter much much longer than 0.005 s - perhaps 0.1 - 0.2 s:

for i1 = 2:2:201,
  subplot(3,1,1)
  plot(t,xPos)                
  hold on
  plot(t,filtfilt(ones(1,i1)/i1,1,xPos),'r')
  hold off
  subplot(3,1,2)
  plot(t,xPos-filtfilt(ones(1,i1)/i1,1,xPos))
  subplot(3,2,5)
  plot(t,gradient(xPos,t))
  subplot(3,2,6)
  plot(t,gradient(filtfilt(ones(1,i1)/i1,1,xPos),t),'r')
  pause(0.1)
end

HTH

[Jan]

+1. You are filtering the positions, while the OP has filtered the velocities. Your approach is much better, because it is focused on the actual noise: the jitter of the determination of the position.

[Marion Rosello]

Thanks for the inputs, Bjorn. The animal actually moves thorughout this acquisition. So, it is (displacement + noise) all along. Coming back to my original point on dampening of the amplitude, no filter size seems to overcome this. I would appreciate if you have any other suggestions.

@Jan - sorry if I was not clear - I was also filtering the displacements in my example above. Only visulazing the velocities later.

[Bjorn Gustavsson]

Yes, I understand that it actually moves, the question is how much of the small-scale stuff is fly-motion and how much is noise. The first thing you have to do is to determine your noise-characteristics, there is no way around that. First you have to look at your data. Really. Just plot the x-position and yposition as a function of time:

subplot(211),plot(t,xPos,'.')
subplot(212),plot(t,yPos,'.')

Zoom in on the first 10 s - clearly your fly is sitting perfectly still during that time-period - both the x and y position gitters between 184.5 and 186 and 45 and 46 respectively with some added smaller noise ontop of that. What you have is discretization noise during that time-period. In your example plot above you have terribly noisy velocities during that first period.

Let's get to the basics... In your example you try to calculate a sensible speed of the fly by filtering the velocities. That is not going to work that well. You have position in 2 dimensions that vary rather slowly except at a couple of instances, but a terrible noise-level. So look at it this way:

A = 1;
B = A/5;
t = (1:5000)/700;
w = 2*pi/3;
xPos = A*sin(w*t) + B*randn(size(t));
dxPos = A*w*cos(w*t);

When you do numerical differentiation of your fly-positions (in this toy example) you get a noise component with approximately this characteristics:

That is your problem. You'll have to estimate a smooth curve whos values can be acceptable as expected values from which your measurements have been obtained through a random process with your excpected noise-characteristics. To obtain such a curve I'd might try something like sgolayfilt over subsets of your time-series with decreasing filterlengths (filter-lengths starting on the order of ~500-1000 samples) and then as soon as the distribution of the residual between the original time-series and the filtered agrees with what you can expect over a a subset, you settle for that filtering over that subset (or even better if it agrees over both the first and second half of that subset.) then continue with increasingly smaller subsets where necessary - this would correspond to a 1-d version of quad-tree decomposition of the signals. Another alternative would be to do the same but with least-square-fitting splines if you have access to the spline toolbox. You can also take a look at John d'Erricos slm-toolbox on the file exchange, there might be similar fitting-tools there you can tailor for this task. Then when you've done the filtering, please make a movie of your fly with the original position and filtered pozsitino plotted with red and blue points - zoomed in to something small enough that we can cearly see the points, say fly-size +/- 10 pixels...

HTH

[Marion Rosello]

Thanks! This helped me think and play more with the data. I'm getting there with the perfect data. :)

[Bjorn Gustavsson]

That's the right attitude! Throw your figurative kitchen sink at the problem, see what gets you forward...

...another thing to try might be to do something in line with what Loren Shure presented about loess smoothing/filtering in one of her post on the mathworks blogs a couple of years ago.

Please let us know how you succeed - because this was an intriguing filtering problem.

All the best