ezfilt

Uses multiple regression based on a selection of cosines.

function [slow, fast] = ezfilt(signal, Fs, cutFreq)

% function [slow, fast] = ezfilt(signal, Fs, cutFreq)

fcrit = cutFreq;

% Create X
X = [];
t = 1:length(signal);
t = (t - 0) ./ Fs;
T = t(end);
f0 = 0;
while 1 == 1,
    cos0 = cos(2 * pi * f0 * t);
    X = [X cos0(:)];
    f0 = f0 + 0.5 / T;
    if f0 > fcrit,
        break;
    end;
end;

% Fit
b = inv(X' * X) * X' * signal(:);
slow = X * b;
fast = signal(:) - slow;

If you have the signal processing toolbox, the normal approach would be to use a real filter as follows:
function vecfilt = teg_filt(vec, Fs, freqLims, hilostr)

% vecfilt = teg_filt(vec, Fs, freqLims, hilostr)
%
% http://www.electronics-tutorials.ws/filter/filter_8.html
% nth order butterworth: H(jw) = 1 / sqrt(1 + eps^2 * (w / w_p)^(2n)))

bw_n = 3; % 6 * bw_n per doubling of frequency
Wn = freqLims ./ (0.5 * Fs);

[z,p,k] = butter(bw_n, Wn, hilostr);
sos = zp2sos(z,p,k);

vecfilt = sosfilt(sos, vec);
vecfilt = sosfilt(sos, vecfilt(end:(-1):1));
vecfilt = vecfilt(end:(-1):1);