I tried using the 'filtic' function (makes initial conditions for 'filter' function) as follows:
% hpf_n = high-pass filter order
% hpf_zi = high-pass filter initial state
% nvar = number of variables
hpf_zi = zeros(hpf_n, nvar);
for i = 1:nvar
hpf_zi(:, i) = filtic(b, a, data(1:hpf_n, i));
end
But it did not initialize the state correctly. The only thing that works is filtering the data backwards and outputting the final state, but this seems unnecessary:
[~, hpf_zi] = filter(b, a, flipud(data));