function fit_spline(obj, f_cut)
	% fit smoothing spline to tracks in obj.tracks
	% evaluate velocity and acceleration in tracks.re_vel and .re_acc
 	% save smoothed position in tracks.re_pred

	n_tracks = numel(obj.tracks);
	% cutoff frequency, rad/s
	omega	= 2*pi*f_cut;
	% spline smoothing parameter
	lambda	= 1/omega^6;

	%% add filtered length field
	b_has_length= any(~isfinite(obj.tracks(1).re_length),'all');
	if b_has_length
		null={nan(3,1,obj.n_t)};
		null=repmat(null, [1 n_tracks]);
		[obj.tracks.fl_length]=null{:};
        [obj.tracks.fl_ldot  ]=null{:};
	end

    %% iterate over all tracks
	for n = 1 : n_tracks
		%% select range of data to interpolate
		% since smoothing spline can only interpolte
        trk		= obj.tracks(n);
		re_ok	= trk.re_pos(:, ~trk.ext);
		t_ok	= trk.t(1, ~trk.ext);

        i_first = find(~trk.ext, 1, 'first');
		i_last  = find(~trk.ext, 1, 'last');
		rng     = i_first : i_last;
		n_interp= numel(rng);
		
		%% filter with smoothing spline
		[re_pred_interp,re_vel_interp,re_acc_interp] = ssplnd(t_ok, re_ok.', trk.t(1,rng), lambda);
		% assign NaN to extrapolated data
		trk.re_pred = nan(3, 1, obj.n_t);
		trk.re_vel  = nan(3, 1, obj.n_t);
		trk.re_acc  = nan(3, 1, obj.n_t);
		trk.re_pred(:,:,rng) = reshape(re_pred_interp.', [3 1 n_interp]);
		trk.re_vel(:,:,rng)  = reshape(re_vel_interp.' , [3 1 n_interp]);
		trk.re_acc(:,:,rng)  = reshape(re_acc_interp.' , [3 1 n_interp]);
		% interpolate data with smoothing spline
		trk.re_pos(:,trk.ext) = trk.re_pred(:,trk.ext);

		%% optional: smooth rod length parameter
		if b_has_length
			re_length	  = trk.re_length(:,~trk.ext);
			[fl_length, fl_ldot] = ssplnd(t_ok, re_length.', trk.t(1,rng), lambda);
			trk.fl_length(:,:,rng) = reshape(fl_length.', [3 1 n_interp]);
            trk.fl_ldot(:,:,rng)   = reshape(fl_ldot.', [3 1 n_interp]);
		end
		
		trk.f_cutoff  = f_cut;

		obj.tracks(n) = trk;
	end
end