filtfilt fails with B=1, A=1 on matrix

filtfilt(1,1,ones(4000,1)); % works OK
filter(1,1,ones(4000,2)); % works OK
filtfilt(1,1,ones(4000,2)); % worked in R2019b Update 6, fails in 2021a Update 5

Error using filter
Initial conditions must be a vector of length max(length(a),length(b))-1, or an array with the leading dimension of size max(length(a),length(b))-1 and with remaining dimensions matching those of x.

Error in filtfilt>ffMultiChan (line 382)
[~,zo] = filter(b(:,ii),a(:,ii),xt,zi(:,ii)*xt(1,:)); % outer product

Error in filtfilt>efiltfilt (line 140)
yCol = ffMultiChan(b2,a2,xCol,zi,nfact,L);

Error in filtfilt (line 89)
y=efiltfilt(b,a,x);

Answer 1

[Walter Roberson]

It appears to be a bug.

In context: filtfilt() calls filter() to do some of the work.

When a and b are scalars, filter is going to be called as

filter(a, b, xt, EMPTY )

where xt is going to have one row and number of columns equal to the number of columns in x. xt is initial conditions of the filters, and the initial values are calculated from the first two rows of the input x.

In this case of scalar a and scalar b, EMPTY is calculated as zeros(0,1) * xt . With xt being size 1 x ncol, the result of the * operation is going to be zeros(0, ncol)

Now...

• when you call filter() with scalar a, scalar b, array X, and EMPTY, then if if EMPTY is 0 x ncol, then filter will succeed for any 2D array size X that is not a row vector (and has ncol columns)
• when you call filter() with scalar a, scalar b, array X, and EMPTY, then if EMPTY is 0 x 1, then filter will succeed for any 2D array size X including row vector X

If you read this carefully, you will see that if X were not a row vector, then having a 0 x ncols initial conditions works. The error message is telling you that your initial conditions must be this size, 0 x ncols. But the error message is wrong for the case where X is a row vector, in which case the initial conditions must be 0 x 1 not 0 x ncols.

If I were to guess, I would guess that X as a row vector is being automatically transposed to X as a column vector, size ncols x 1, and it is wanting the 0 x ncols to match that number of columns (1) .

If I am correct, then the problem would be fixed if filtfilt() were to pass in the dimension to filter along.

In the meantime, all you can do is filtfilt() one column at a time of your signal.

[Alain]

Thanks for the explanation. I'm not sure I follow all the details. X was a column vector.

More generally, I'd expect (a) code that worked in one version not to fail in the next, (b) code that works in the general case (B and or A vectors) to work in a limit case (B and A scalars), and (c) not to have to rewrite code to work around the problem.

Setting B=1; A=1; is a simple and clean way to switch off filtering, easier than adding tests to make filtering conditional on a flag. Rewriting code is (a) time consuming, (b) error prone, (c) impractical for code that has shipped, and (d) no guarantee that the issue is really solved given that it arose for mysterious reasons. We really need to be able to trust the basics.

[Walter Roberson]

Your ones(4000,2) is not a vector . If it had been a vector, filtfilt() would have detected it and used a different internal routine.

In the routine that is used for the case where multiple channels (columns) exist, because your filter lengths are both 1, the current version of filtfilt() ends up passing a row vector to filter(), with the intention that each column be filtered separately. However, as is common for a lot of functions, when filter() receives a row vector, it assumes that it represents a single signal, instead of assuming that it is receiving multiple channels of a signal that has only one sample per channel.

If your filter length had been larger than 1, then something with multiple rows would have been passed in to filter, and everything would have been fine.

You accidentally triggered an edge case. Now that I have sent a bug report with analysis to Mathworks it should be simple for Mathworks to fix the code by passing in the dimension to operate along (1, that is)

More generally, I'd expect (a) code that worked in one version not to fail in the next,

Yes, understandable, valid code is not expected to fail, not unless the release notes document a change

(b) code that works in the general case (B and or A vectors) to work in a limit case (B and A scalars

No, this is not an valid assumption. MATLAB has a lot of places where the code works along the first non-scalar dimension unless you tell it which dimension to work along. Unfortunately, filtfilt() does not offer an option to pass in the dimension to work along.

and (c) not to have to rewrite code to work around the problem.

Ummm -- I am not sure if you are aware, but I do not work for Mathworks. I am the volunteer who stayed up all night and debugged your code at about 5 AM my time.

By the way: you can set B = [1; 0] A = 1 to turn off filtering with multi-channel inputs.

[Alain]

Sorry, I didn't intend to criticize you. My point was that the issue should be fixed (by MW if listening), rather than just explained and/or worked around. Thanks again for your input!

[Walter Roberson]

Mathworks mostly doesn't read these postings. Some individual Mathworks employees volunteer to answer questions that interest them... if they happen to notice the posting.

The way to get things fixed is to open a Support case about it. Which I have done for this already, and pointed them to this topic for backing information.

If you use the Contact form on my profile to send me your email address, I will ask Mathworks to add you to the case to be notified when it is fixed.