Getting NaN after filtfilt function

I am applying a simple butterworth bandpass filter to a signal having low cut off freq 0.5 and higher cutoff freq 3. But after designing the filter when i am applying filtfilt function i am getting the signal as array of NaN. To check if in the specified frequency range any signal does not exist, i crosschecked it with fft and there was some peaks in that range, which indicate the presence of signl in that range. I dont understand why i am getting NaN.

I am giving the data link of matlab drive, fft plots and code.

Matlab Drive Link of DATA

sampling_frequency = 1000000;
signal = concatenatedData(:, 2); 
order = 11; 
low_cutoff = 0.5; 
high_cutoff = 3;  
[b, a] = butter(order, [low_cutoff, high_cutoff] / (sampling_frequency / 2)); 
filtered_signal = filtfilt(b, a, signal); 

Answer 1

[Star Strider]

You are using transfer function (‘[b,a]’) output implementation. Often, this creates an unstable filter, producing NaN values in the fitlered output. Using freqz, the instability is readily apparent —

sampling_frequency = 1000000;
% signal = concatenatedData(:, 2); 
order = 11; 
low_cutoff = 0.5; 
high_cutoff = 3;  
[b, a] = butter(order, [low_cutoff, high_cutoff] / (sampling_frequency / 2)); 
% filtered_signal = filtfilt(b, a, signal); 

figure
freqz(b, a, 2^16, sampling_frequency)
set(subplot(2,1,1), 'XScale','log')
set(subplot(2,1,2), 'XScale','log')

The solution is to use second-order-section implementation instead, as discussed in the butter documentation section on Limitations. (This is true generally, not simply with Butterworth filters, and the same discussion appears in the documentation on all of them.)

sampling_frequency = 1000000;
% signal = concatenatedData(:, 2); 
order = 11; 
low_cutoff = 0.5; 
high_cutoff = 3;  
[z, p, k] = butter(order, [low_cutoff, high_cutoff] / (sampling_frequency / 2)); 
[sos,g] = zp2sos(z,p,k);
% filtered_signal = filtfilt(sos, g, signal); 

figure
freqz(sos, 2^16, sampling_frequency)
set(subplot(2,1,1), 'XScale','log')
set(subplot(2,1,2), 'XScale','log')

There is no instability in the second-order-section implementation.

EDIT — (15 Jan 2024 at 16:52)

Added the buttord call, that designs a 113-order filter —

Wp = [low_cutoff, high_cutoff] / (sampling_frequency / 2);
Ws = [0.9 1.1] .* Wp;
Rp =  1;
Rs = 60;
[n,Wn] = buttord(Wp,Ws,Rp,Rs)

n = 59

Wn = 1×2

1.0e-05 * 0.0991 0.6053

[z,p,k] = butter(n,Wn);
[sos,g] = zp2sos(z,p,k);
% filtered_signal = filtfilt(sos, g, signal); 

figure
freqz(sos, 2^16, sampling_frequency)
set(subplot(2,1,1), 'XScale','log')
xt1 = get(subplot(2,1,1), 'XTick')

xt1 = 1×3

1000 10000 100000

set(subplot(2,1,2), 'XScale','log')
xt2 = get(subplot(2,1,2), 'XTick')

xt2 = 1×5

10 100 1000 10000 100000

NOTE — A passband of [0.5 3.0] Hz with a 1 MHz sampling frequency may be difficult to realise in the best of circumstances.

.

[Ritesh]

Thank you so much for the information. As you mentioned in the note, A passband of [0.5 3.0] Hz with a 1 MHz sampling frequency may be difficult to realise in the best of circumstances. I am curious if there is any method to know the value of sampling frequency and corresponding minimum frequency range to realize a bandpass successfully.

[Star Strider]

As always, my pleasure!

That is most likely a problem of frequency resolution in both the time domain and frequency domain. A 3 Hz passband with at 500 KHz Nyquist frequency equates to as the normalised frequency. This is going to be difficult to realise in almost any filter design. Different filter designs may be more successful than others, for example an elliptic design can manage that (illustrated here), however the filtered result might not be acceptable, since filtering it would eliminate too much detail in the signal. Note here that I needed to specify a Fourier transform length of 2^22 = 4194304 to get this level of detail with freqz. This code designs a 6-order elliptic filter, the 133-order Butterworth design was not capable of this sort of detail when I ran it using similalr code.

Unless you need a maximally flat passband and a much more gradual stopband, consider using this elliptic filter in your current application —

Fs = 1E+6;
Fn = Fs/2;
Wp = [0.50 3.0]/Fn;                                                 % Passband Frequency (Normalised)
Ws = [0.10 4.0]/Fn;                                                 % Stopband Frequency (Normalised)
Rp =  1;                                                            % Passband Ripple
Rs = 60;                                                            % Passband Ripple (Attenuation)
[n,Wp] = ellipord(Wp,Ws,Rp,Rs);                                     % Elliptic Order Calculation
[z,p,k] = ellip(n,Rp,Rs,Wp);                                        % Elliptic Filter Design: Zero-Pole-Gain 
[sos,g] = zp2sos(z,p,k);                                            % Second-Order Section For Stability

figure
freqz(sos, 2^22, Fs)                                                % Filter Bode Plot
set(subplot(2,1,1), 'XScale','log', 'XLim',[0 10])
set(subplot(2,1,2), 'XScale','log', 'XLim',[0 10])

Experiment with the upper stopband frequency (‘Wp’) to get the result you need. (I would not change the lower stopband frequency, since it is too close to zero to be extended, and increasing it much may not produce a desirable filtered result, since the transition region would be too narrow.) Changing the stopband attenuation (‘Rs’) is another option, again, decreasing it would likely produce a better filtered result than increasing it.

Your design choices depend on what you want to do and the filtered result you are willing to accept.

.

[Ritesh]

Thak you. I got the point.

[Star Strider]

As always, my pleasure!

[Ritesh]

Hello Star Strider. As it is generally hard to implement a BPF of such low frequency bands, would it be a good idea to downsample my original signal to a suitable frequency.

[Paul]

Hi Ritesh,

Taking a look at your original filter, with @Star Strider suggestion to implement via sos.

sampling_frequency = 1000000;
order = 11; 
low_cutoff = 0.5; 
high_cutoff = 3;  
[z, p, k] = butter(order, [low_cutoff, high_cutoff] / (sampling_frequency / 2)); 
[sos,g] = zp2sos(z,p,k);

With such a high sampling frequency, the number of points at which to compute the frequency response needs to be much higher than 2^16. Using 2^16 points, the resolutio of freqz will be

sampling_frequency/2/2^16  % Hz

ans = 7.6294

The first frequency point would be too large to see the action below 5 Hz. Use many more points, so that we get 0.5 Hz resolution

%freqz(sos, 2^16, sampling_frequency)
[h,f]=freqz(sos, sampling_frequency, sampling_frequency);

figure
plot(f,db(abs(h*g))) % multiply by g to get correct gain

Zoom in to see what's happening at low frequency, don't plot in dB to see what's happening at f = 0.

figure
plot(f,(abs(h*g)),'-o'),xlim([0 20])
xline([low_cutoff, high_cutoff])

If we really want to see what's going at around the passband, use a lot of points in the call to freqz just around the pass band.

[h,f]=freqz(sos, 0:.001:10, sampling_frequency);

figure
plot(f,(abs(h*g))) % multiply by g to get correct gain
xline([low_cutoff, high_cutoff])

[Star Strider]

@Ritesh — The problem is the frequency resolution, and actually the time resolution as well (since the filter will actually be operating in the time domain). I am not certain that any filter could work effectively in this sort of application (considering the difference between the passband frequencies and the sampling frequency), however the much more computationally efficient elliptic design can do much better here than a Butterworth design, as my example demonstrates. Unless you absolutely need a maximally flat passband, an elliptic filter is an appropriate choice. The passband ripple in the filter I demonstrated is only 1 dB, so it approaches a flat passband.

As for downsampling to a lower sampling frequency, that is a tradeoff you will need to consider. Downsampling will by definition eliminate some signal details that would be present in the original signal, and if those details were present in the filter passband, they would be lost. If they were preserved in the passband, downsampling might be appropriate.

One way to determine that would be to calculate the Fourier transforms of the original signal and of the downsampled signal, and look at the signal characteristics in the filter passband range of 0 to 5 Hz. Choose a Fourier transform length that is identical for both of them (the second argument in the fft function), so that you can accurately compare the two spectra, since in that event they will have the same frequency resolution. If they are essentially identical, then nothing would be lost to downsampling. If they are significantly different (and only you can determine that), downsampling might not be appropriate.

This is essentially a theoretical discussion, since I have not encountered a situation similar to this (none of the signals I work with are sampled at greater than about 50 kHz), so I have no direct experience with it. I am not even certain how to simulate it in order to test these ideas. The best I can do is to outline an approach to determine if downsampling would be appropriate.

[Paul]

I'm confused as to whether or not there really is a problem. Why would the difference between the passband frequencies and the sampling frequency be an issue? As shown above, the frequency response of the filter seems to be as expected, at least insofar as the width of the passband. Here's an example of filtering a simple signal in the time domain that appears to give the expected result.

Define the filter

sampling_frequency = 1000000;
order = 11; 
low_cutoff = 0.5; 
high_cutoff = 3;  
[z, p, k] = butter(order, [low_cutoff, high_cutoff] / (sampling_frequency / 2)); 
[sos,g] = zp2sos(z,p,k);

Define the signal as the sum of two sine waves at 2 and 10 Hz respectively

F1 = 2;
F2 = 10;
t = 0:(1/sampling_frequency):10;
u = sin(2*pi*F1*t) + sin(2*pi*F2*t);

Filter the input. I arbitrarily chose to mulitply the ouput of sosfilt rather than the input to sosfilt.

y = g*sosfilt(sos,u);

Plot the result and compare to expected steady state output based on the frequency response.

h = freqz(sos, [1 2 3], sampling_frequency);
phi = angle(h(2));
figure
plot(t,y,t,sin(2*pi*F1*t+phi)),grid

As should be the case, essentially only the 2 Hz component shows up in the steady state output.

One thing I did notice is that the value of g is quite small

g

g = 7.0140e-57

I guess that means that the peak output of sosfilt is quite large, i.e., on the order of 1/g. I don't know if there are any implications of that on numerical accuracy of sosfilt, or if it would be better to multiply the input by g, before applying sosfilt, or if some other scaling is in order. I do see that zp2sos offers additional inputs to control the ordering and coefficient scaling of the individual sos sections, but I don't if those additional inputs would be helpful in this case.

[Ritesh]

Hello @Paul. you mean if my frequency of interest for BPF is too small as compared to the sampling frequency, there should not be any issue? But I was getting NaN values for filtered signals. At least it can not be handled in a direct way. In the previous example you have taken a sinusoid of 2 and 10Hz, but in my case the input signal is sampled at 1000000Hz and it contains many frequency components. I am trying to extract the signal in the range of 0.5 to 3Hz.

[Paul]

As has already been stated, you were likely getting NaN values for your filtered signals because you were using the transfer function form (b,a) of the filter rather than the sos form. Even though the tf and sos forms are the same mathematically, the latter is supposed to be much more robust numerically, particularly for higher order filters. This idea is probably discussed somewhere in the documentation.

Instead of providing a link to Matlab Drive, why not upload an example input signal in a .mat file (use the paperclip icon on the Insert menu) and copy/paste to here the code that illustrates the NaN problem.

I'm not sure why having your passband at a very low frequency compared to the sampling frquency should, in and of itself, be an issue, at least mathematically. Yes, I showed the simple sum of two sinewaves, but they were both sampled at the sampling_frequency you provided and the sos filtering seemed to work fine. Post your input signal(s) and we can see if the sos filtering works for you.

[Ritesh]

Hello @Paul. Sorry for late reply. Actually my file file is in GB so it is difficult to upload. The data is of radar. I wanted to extract a particular frequency corresponding to a range, to check if it is working. So what I found that it is better to plot the fft in a frequency range instead of extracting the signal in time domain. While i am seeting xlim in sich small range plot is not that much good due to resolution issues. So I have increased the frequency resolution by dividing the sampling frequency with 10000*N instead of only N. Now I am getting the expected peak. I was wondering if it is the correct method.

signal = signal - mean(signal);
sampling_frequency = 1000000; 
N = length(signal); 
fft_result = fft(signal); 
frequencies = (0:N-1) * (sampling_frequency/(10000*N)); % Frequency values
figure; 
plot(frequencies, abs(fft_result)); 
title('FFT'); 
xlim([0.5 3])
xlabel('Frequency (Hz)');  
ylabel('Amplitude'); 

[Star Strider]

You can easily downsample a signal using the resample function. The only problem is that you lose the high-frequency details in the signal. If there are no high-frequency details that you want or need, then downsampling presents no problems. The filter must, of course, be re-designed to the new sampling frequency.

EDIT — (27 Jan 2024 at 15:02)

I was not certain how downsampling might afffect the signal (downsamplilng is a form of lowpass filtering, since it essentially removes frequencies above the downsampling frequency) so I decided to explore downsampling and the highpass filtering.

These are the results —

Fs = 1E+6;
L = 25;
t = linspace(0, Fs*L, Fs*L+1)/Fs;
Fs_check = 1/(t(3)-t(2))

Fs_check = 1000000

fm = 0.25:0.25:8;
s = sum(sin(2*pi*fm(:)*t));

[FTs1,Fv] = FFT1(s,t);
                                                                                                                                                                                                    figure
plot(Fv, abs(FTs1)*2)
grid
xlim([0 10])
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('Original')

Fs = 1E+6;
Fn = Fs/2;
Wp = [0.50 3.0]/Fn;                                                 % Passband Frequency (Normalised)
Ws = [0.10 4.0]/Fn;                                                 % Stopband Frequency (Normalised)
Rp =  1;                                                            % Passband Ripple
Rs = 60;                                                            % Passband Ripple (Attenuation)
[n,Wp] = ellipord(Wp,Ws,Rp,Rs);                                     % Elliptic Order Calculation
[z,p,k] = ellip(n,Rp,Rs,Wp);                                        % Elliptic Filter Design: Zero-Pole-Gain 
[sos,g] = zp2sos(z,p,k);                                            % Second-Order Section For Stability

figure
freqz(sos, 2^22, Fs)                                                % Filter Bode Plot
set(subplot(2,1,1), 'XScale','log', 'XLim',[0 10])
set(subplot(2,1,2), 'XScale','log', 'XLim',[0 10])

                                                                                                                                                                                                                                                                                                
s_filt = filtfilt(sos,g,s);

[FTs_filt1,Fv] = FFT1(s_filt,t);

figure
plot(Fv, abs(FTs_filt1)*2)
grid
xlim([0 10])
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('Bandpass Filtered')

Fsr = 6;                                                        % NOTE — Two Times The Desired High Passband Frequency
[sr,tr] = resample(s, t, Fsr);

[FTsr1,Fv] = FFT1(sr,tr);

figure
plot(Fv, abs(FTsr1)*2)
grid
xlim([0 10])
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('Downsampled')

sr_filt = highpass(sr, 0.5, Fsr, 'ImpulseResponse','iir');

[FTsr_filt1,Fv] = FFT1(sr_filt,tr);

figure
plot(Fv, abs(FTsr_filt1)*2)
grid
xlim([0 10])
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('Downsampled and Highpass Filtered')

function [FTs1,Fv] = FFT1(s,t)
    s = s(:);
    t = t(:);
    L = numel(t);
    Fs = 1/mean(diff(t));
    Fn = Fs/2;
    NFFT = 2^nextpow2(L);
    FTs = fft((s - mean(s)).*hann(L), NFFT)/sum(hann(L));
    Fv = linspace(0, 1, NFFT/2+1)*Fn;
    Iv = 1:numel(Fv);
    FTs1 = FTs(Iv);
end

I am simply presenting this as a simulation. I am not making any specific claims about it.

.

[Ritesh]

@Star Strider thanks. Got to know new things. Downsampling is nothing but LPF. Signal processing is so interesting. :-) As the discussion is going on. I want to learn more. So how far I can go with downsampling? like if i am only looking for 3 to 5 hz range. Can i donsample it from 1Mhz to 100Hz or so? Or will there be any issue with that.

I know my issue can be handled with resample now. However just curious to know about my previous code "frequencies = (0:N-1) * (sampling_frequency/(10000*N)); ". here if I am increasing the resolution by dividing fs with high number then i am also disturbing the time domain signal, something like zero padding.

[Paul]

What exactly was this line

frequencies = (0:N-1) * (sampling_frequency/(10000*N));

supposed to accomplish?

Those frequencies do not correspond to the output of this line:

fft_result = fft(signal);

Of course, you can make the plot

plot(frequencies, abs(fft_result));

but it's not a correct representation in the frequency domain.

[Star Strider]

@Ritesh —

My pleasure!

I am not certain what you are doing with the ‘frequencies’ assignment. It is posible to downsample by simple indexing, however that approach is not recommended. The resample function uses an anti-aliasing fitler so that unwanted frequencies are filtered out. I would use if tor that reason.

My example code simply demosntrates the effect of resampling and the effect of filtering a resampled signal.

[Ritesh]

Ok thanks for the clarification.

[Star Strider]

As always, my pleasure!

Answer 2

[Hassaan]

1. Stability of the Filter: With high filter orders, like the 11th order you're using, the filter can become unstable, especially when applied in a forward-backward manner as filtfilt does. The Butterworth filter design should theoretically be stable, but numerical issues can arise due to finite precision arithmetic.
2. Initial Conditions: The filtfilt function typically handles initial conditions to reduce transients by using a certain method of padding. If something is off with the initial conditions or the padding, it might lead to NaN values.
3. Data Quality: Ensure there are no NaN or Inf values in your input signal. filtfilt will not be able to handle these and will propagate NaNs throughout the output.
4. Filter Coefficients: Check the filter coefficients b and a for any NaN or Inf values. This could happen if the butter function encounters numerical issues.
5. Signal Characteristics: If the signal has certain characteristics that are not compatible with the filtering process, it might lead to issues. However, since you mentioned that the FFT shows peaks within the passband, this seems less likely.

Here's what you can do to debug:

• Check for NaNs or Infs in your input signal and filter coefficients.
• Reduce the filter order. A lower order might be more stable numerically.
• Normalize your signal if it has very large or very small values. Such values could cause numerical instability.
• Use a different filtering approach or design method to see if it's specific to filtfilt or Butterworth design.

This script includes checks for NaNs in your signal, filter stability, and also reduces the filter order for better numerical stability.

% Load your signal data
% Assuming 'concatenatedData' is already in your workspace
signal = concatenatedData(:, 2); 

% Check for NaNs or Infs in the input signal
if any(isnan(signal)) || any(isinf(signal))
   error('Signal contains NaNs or Infs.');
end

% Define filter parameters
sampling_frequency = 1000000;
order = 2; % Reduced order for stability
low_cutoff = 0.5; 
high_cutoff = 3;  

% Calculate filter coefficients
[b, a] = butter(order, [low_cutoff, high_cutoff] / (sampling_frequency / 2), 'bandpass');

% Check for NaNs or Infs in filter coefficients
if any(isnan(b)) || any(isnan(a)) || any(isinf(b)) || any(isinf(a))
   error('Filter coefficients contain NaNs or Infs.');
end

% Check stability of the filter
if any(abs(roots(a)) >= 1)
   error('Filter is unstable.');
end

% Apply the filter using filtfilt
filtered_signal = filtfilt(b, a, signal);

% Check for NaNs in the filtered signal
if any(isnan(filtered_signal))
   nan_indices = find(isnan(filtered_signal));
   fprintf('NaNs start at index %d.\n', nan_indices(1));
else
   fprintf('Filtering completed without NaNs.\n');
end

% Plot the original and filtered signals for comparison
figure;
subplot(2,1,1);
plot(signal);
title('Original Signal');
subplot(2,1,2);
plot(filtered_signal);
title('Filtered Signal');

---------------------------------------------------------------------------------------------------------------------------------------------------------

If you find the solution helpful and it resolves your issue, it would be greatly appreciated if you could accept the answer. Also, leaving an upvote and a comment are also wonderful ways to provide feedback.

Professional Interests

• Technical Services and Consulting
• Embedded Systems | Firmware Developement | Simulations
• Electrical and Electronics Engineering

Feel free to contact me.

[Ritesh]

Thank you Muhammad for the information. I came to know that the filter is unstable. The reason is clear with star striders answer.