filtfilt function returning NaN at certain frequencies

I'm hoping someone can help me sort this out. I'm trying to use the filtfilt function with a butterworth filter, and the maximum frequency range I can specify depends on the sampling frequency Fs. For example:

Fs=16kHz, my maximum bandpass range is [121 7944]. Any frequency lower than 121 or higher than 7944 causes filtfilt to return NaN.

Fs=32kHz, my maximum bandpass range is [209 15900]. Any frequency lower than 209 or higher than 15900 causes filtfilt to return NaN.

Similar questions have been asked in the forum before, but the answers were usually to check the data for NaN or Inf values. I have neither in my data.

The program I'm developing needs to be generic enough to accept any sampling rate and user-specified bandpass. I want to add a warning if the user specifies a bandpass that will cause problems, but I'm not sure how to calculate the limit. I'm sure a large part of my problem is unfamiliarity with filters. Is there some nuance in filter design that I'm missing? Any help would be appreciated. Thanks.

Michael

[Ryan]

Did you ever find a solution to your problem? I'm running into the very same. I can filter at high frequencies, but my target bandpass is 8-30Hz and all I get are NaNs. Checked and double checked for inf and NaN in the original dataset and nothing.

[Jan]

@Ryan: Are you using the Octave implementation of FILTFILT as included in EEGLab also? Then you could try FILTFILTM from http://www.mathworks.com/matlabcentral/fileexchange/32261-filterm, which is compatible to Matlab's version, but faster when the MEXed replacement of filter is used.

As you see, I cannot post this as an answer, although it is one, because you did not post your question as a question, but as a comment. But I see the reason for your decision. Perhaps you could provide some example data such that the problem can be reproduced.

[Jeremy]

I am seeing this same issue and it is not related to input as I can reproduce it with randn numbers:

x=randn(35000,1);
[b, a] = butter(10, 60/5000, 'high');
x=filtfilt(b,a,x);

[Jan]

@Jeremy: And again, which version of FILTFILT are you using? Matlab < or >= 2011b, Octave's version, the implementation in EEGLAB or the compiled MEX from the FileExchange? I couldn't reproduce any of the problems in this thread because in every description a single important detail is missing. But I'm really eager to solve this.

[Jeremy]

2012b, 64-bit. I am just using the built-in butter and filtfilt functions.

edit: My neighbor ran this on 2012b 32-bit and got a valid result, I ran on 2010b 64-bit and still got an invalid result (not NaNs, but numbers in the -10^60 to +10^144 range.

Answer 1

[Cameron]

I found the source of my problem. I am using the EEGLab toolbox. It appears that it causing a problem with the native filtfilt fxn. Removing the EEGlab toolbox from my paths list fixed the problem.

Answer 2

[Jan]

Reply to Jeremy's comment:

x = randn(35000,1);
[b, a] = butter(10, 60/5000, 'high');
x = filtfilt(b,a,x);

Inside filtfilt the reflection method for the start values need to solve the linear system (from version 2009a):

zi = sp \ (b(2:nfilt).' - a(2:nfilt).'*b(1));

But for the filter parameters a and b of the Butterworth filter, the condition of sp is 2.3e19. As a rule of thumb you loose k significant digits in the solution of the linear system, when the matrix has the condition 10^k. In this case loosing 19 digits mean, that the results are random. Huge values and overflows can be expected.

Obviously the reflection method proposed by Fredrik Gustafsson (Determining the initial states in forward-backward filtering, IEEE Transactions on Signal Processing, pp. 988--992, April 1996, Volume 44, Issue 4) has a limitation here. Unfortunately I do not have access to this paper to find out, what is said in the discussion about this problem. Sending an enhancement request to TMW would be a good idea.

The FiltFiltM implemtation suffers from exactly the same problem, but it does not use a sparse method to solve this system, such that a warning message appears. I'm not sure, why this does not happen for the sparse version.

[Jeremy]

Thanks! any idea why 64-bit Matlab is more susceptible? I can cause it on 32-bit Matlab but it takes a cutoff frequency closer to the extremes.

My problem is that I am implementing this in a validated GUI, so if there are limits I need to know what those are so I can error trap those conditions. So far there doesn't seem to be any real definite limit. I have submitted a service request so I'll see where that goes.

[Jan]

You could use a modified copy of FILTFILT or FiltFiltM from the FileExchange, which check the condition of the reflection matrix at first. When the condition is above 10^10, the method to reduce transitional effects at the start and end of the signal will fail. A simple solution is to disable the tricks to reduce the transitional effects. A warning should appear at least, but it is not a bug to admit, that the the filter parameters do not allow to apply trials to suppress noise induced by the filter algorithm. Physical filters e.g. in electric circuits suffer from the same effects.

For FiltFiltM:

K = ...
if condest(K) < 1e10
  IC = K \ (b(2:Order) - a(2:Order) * b(1));
else
  IC = 0;
  warning('JSimon:FiltFiltM:TransitionalNoise', ...
          'Cannot damp transitional effects at the edges.'); 
end

[Jeremy]

Talked to TMW, the workaround is to use the filter spec object, seems to work well. No issues with long records or extreme cutoff frequencies.

d = fdesign.highpass('N,F3dB',10,60/5000);
hd = design(d,'butter');
x=randn(35000,1);
y = filtfilt(hd.sosMatrix,hd.ScaleValues,x);

Answer 3

[Jan]

filtfilt is based on the calculations of filter, which is a deterministic function. When you call it with finite values, infinite values are only replied in case of overflows.

I would be helpful if you post the relevant part of the code instead of describing it by some text. E.g. your maximum bandpass range of [121 7944] looks mysterious. How do you define the filter parameters exactly? Usually the cut-of frequencies are normalized by the Nyquist frequency, such that only values between 0 and 1 are accepted. Perhaps your test for finite values of the input has a bug and inspite of your assumption the input does contain NaNs.

Anyhow, whenever I've seen reports of such magic behaviour, there have been very real bugs at unexpected locations. Usually Matlab's toolbox functions have been shadowed by user-defined functions or variables. And without any doubt, Matlab's filter accepts the full range of possible frequencies from 0 to 1 (exclusively for numerical reasons). You can see an implementation as M-file with exactly the same results here: Answers: Use Filter Constants. So you can use this M-function instead of Matlab's built-in filter() and check in every step, where any why the non-finite values occur.

[Michael]

Thanks so much for the quick response Jan. Sorry - I wasn't very clear in the original post. My code is:

lowerFreq=handles.mfaPrefs.fmin*2/handles.Fs;
upperFreq=handles.mfaPrefs.fmax*2/handles.Fs;
[b a]= butter (10,[lowerFreq upperFreq]);
xxx=filtfilt(b,a,x);

where handles.Fs is the sampling frequency, handles.mfaPrefs.fmin is the lower frequency limit, handles.mfaPrefs.fmax is the upper frequency limit, and x contains the data.

I came up with the values in the original post through trial and error. When Fs=16000, I could only set fmin=121 and fmax=7944. That puts the vector passed to the butter function at [.0151 .993]. If I lowered fmin to 120 or increased fmax to 7945, the filtfilt function returned a vector of NaN.

I don't have any NaN or Inf values in the dataset. Both find(isnan(x)) and find(isinf(x)) return empty matrices The data ranges from -4.88E-4 to .0278.

I just tried it, and the filter() function seems to have no problem with this. I changed the range to [.001 .999] and still got real numbers. But I'm worried about phase shifts using the filter function. I'm analyzing acoustical data, and logging the time of certain acoustic events. Phase shifts would cause problems.

Any suggestions would be appreciated. I'm a little stuck.

Michael

[Michael]

One quick follow-up. I stepped through the filtfilt function and the problem lies with this line:

y = [2*x(1)-x((nfact+1):-1:2);x;2*x(len)-x((len-1):-1:len-nfact)];

According to the comments, this command is used to Extrapolate beginning and end of data sequence using a "reflection method". Slopes of original and extrapolated sequences match at the end points. This reduces end effects.

If I run my original dataset (x) through the rest of the code, no problem. But when it's transformed into y, the first filter() call returns NaN for over half the dataset.

[Jan]

@Michael: This sounds such strange, that I believe there is another bug anywhere else.

Answer 4

[Cameron]

I am having the same problem. I am using OSX 10.8.2 with 2012a. filtfilt was working fine but now using the 'butter' filter is only returning NaN's. If I used the code in the example found in doc filtfilt the 'equiripple' works fine but the the 'butter' again returns garbage. I wonder is this has something to do with a Java update or OSX Update.

doc filtfilt

rng default;
x=ecg(500)'+0.25*randn(500,1); %noisy waveform
h=fdesign.lowpass('Fp,Fst,Ap,Ast',0.15,0.2,1,60);
d=design(h,'equiripple'); %Lowpass FIR filter
y=filtfilt(d.Numerator,1,x); %zero-phase filtering
y1=filter(d.Numerator,1,x); %conventional filtering
subplot(211);
plot([y y1]);
title('Filtered Waveforms');
legend('Zero-phase Filtering','Conventional Filtering');
subplot(212);
plot(ecg(500));
title('Original Waveform');

h=fdesign.lowpass('N,F3dB',12,0.15);
d1 = design(h,'butter');
y = filtfilt(d1.sosMatrix,d1.ScaleValues,x);
plot(x,'b-.'); hold on;
plot(y,'r','linewidth',3);
legend('Noisy ECG','Zero-phase Filtering','location','NorthEast');

[Malcolm Lidierth]

Works fine for me with:

Operating System: Mac OS X Version: 10.8.2 Build: 12C60 Java Version: Java 1.6.0_41-b02-445-11M4107 with Apple Inc. Java HotSpot™ 64-Bit Server VM mixed mode --------------------------------------------------------------------------------------------------------

[Jan]

Does dbstop if NaNInf reveal the line, which creates the NaNs?

Answer 5

[Cameron]

Right, that is my point. It used to work for me and then post this last OSX update where apple mucked around with Java because of the security issue it is no longer working.

MatLab Version 7.14.0.739 (R2012a) OSX 10.8.2

java version "1.7.0_05" Java™ SE Runtime Environment (build 1.7.0_05-b06) Java HotSpot™ 64-Bit Server VM (build 23.1-b03, mixed mode)

Answer 6

[Carlo E.D. Riboldi]

I had a very similar issue on dataset sampled at very high frequency and very long. NaN showed up without apparent regularity when applying filtfilt to the signal, changing the order of the Butterworth filter (butter) between 4 and 8 and the cut-off frequency.

I tried the implementation of the filter through filter objects (fdesign), but this didn't fix the problem.

I downloaded the small package FilterM from Matlab Central, and this apparently fixed it all (<http://www.mathworks.com/matlabcentral/fileexchange/32261-filterm)>. I recompiled the included Mex-file with MSVC2008.

I am using Matlab 7.9.0, R2009b.