I'd say this definitely needs to be looked at.  On my system I had a filtfilt.m version in my path from probably 5+ years ago and it is extremely simple, i.e., essentially calling filter in both directions after tagging some zeros at the end of the input signal to capture the filter transient and use that as the start of the reverse filter.  Here is the comment from that old version:

## TODO: In Matlab filtfilt `reduces filter startup transients by carefully
## TODO:    choosing initial conditions, and by prepending onto the input
## TODO:    sequence a short, reflected piece of the input sequence'.
## TODO:    Once filtic is written, use that here.

There is a filtic() function, but I don't see it called anywhere in the latest filtfilt.m.  Not sure it applies.

I experimented a bit, and I found that the main contributing factor is this line:

    v = [2*x(1,c)-x((lrefl+1):-1:2,c); x(:,c);
         2*x(end,c)-x((end-1):-1:end-lrefl,c)]; # a column vector

If I change that line to the following:

    v = [zeros(size(2*x(1,c)-x((lrefl+1):-1:2,c))); x(:,c);
         zeros(size(2*x(end,c)-x((end-1):-1:end-lrefl,c)))]; # a column vector

the transient at the front end goes away.  I guess that is no surprise, but why the 2 of 2*x(1,c)?  Consider if filter length is 2

x(1,1) = 3
x(2,1) = 1.5
x(3,1) = -0.5
...

Then

v(1) = 3 + 3 - (-0.5) = 6.5
v(2) = 3 + 3 - (1.5)  = 4.5
v(3) = 3
v(4) = 1.5
v(5) = -0.5

The point is that there is always some additional bias in this padding meant to prime the filter.

I'm not following the logic of using si*v(1) for the initial state either.  I took out the padding at the front and still see a small transient that isn't at the level of the random signal.  Guess I'd have to read the paper referenced within, but I don't see how this initialization scheme makes sense.

Further note, your example runs very slow because filtfilt.m is using a for-loop whereas filter might not be.  It is clear that the F = flipud ( filter ( b, a, flipud ( filter ( b, a, M ) ) ) ); line is faster than filtfilt.

Filtfilt may introduce strange artefacts into the spectra and time series of filtered signals.

This seems to result from the code computing the initial filter state, and manifests itself as a ringing artefact at both ends of the filtered signal, which does not appear when simply applying "filter" in both directions.

See the attached file for sample code plotting a short demo.

While being far from catastrophic for most applications, the effect is certainly strange enough to cause a certain amount of headache to unprepeared users, as it does not seem to be documented anywhere, and it might break assumptions regarding expectation values towards both ends of time series.

Is this really how we want in this function to work?