Filtfilt returns NaN matrix

Question

Aboltabol 2023-5-10

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1961079-filtfilt-returns-nan-matrix

评论： Star Strider 2023-5-11

I am running the following operation (please download the attached file to your working directory):

load('Channel_Sim.mat');
[b,a] = butter(4, [3 20] ./ (1000/2)); % Sampling frequency is 1000 Hz.
chan_filtered = filtfilt(b, a, chan_data)

However, chan_filtered is a NaN matrix. Why?

I have checked that chan_data does not contain any NaN or Inf. I also tried resetting the butterworth filter range (3-20) over a wide range of values but to no avail. A PSD plot (see attached jpeg file) shows that chan_data encompasses a wide range of frequencies inluding the target frequency (3-20).

2 个评论
显示无隐藏无

Walter Roberson 2023-5-10

If you have even 1 nan or inf in your data the filtered results will likely be nan

Aboltabol 2023-5-10

在 MATLAB Online 中打开

@Walter Roberson : I have checked that.

>> sum(isnan(chan_data))
ans =
     0
>> sum(isinf(chan_data))
ans =
     0

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2023-5-10

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1961079-filtfilt-returns-nan-matrix#answer_1232589

在 MATLAB Online 中打开

Channel_Sim.mat

The transfer function implementation of the filter is unstable, however I cannot determine the reasson for the instability, although it may simply be the relative magnitudes of the numerator and denominator coefficients.

The easiest way to solve that problem is to use second-order-section implementation of the same filter —

load('Channel_Sim.mat')

% chan_data

Fs = 1000;

[b,a] = butter(4, [3 20] ./ (1000/2));

[z,p,k] = butter(4, [3 20] ./ (1000/2));

[sos,g] = zp2sos(z,p,k);

sos = 4×6

1.0000 2.0000 1.0000 1.0000 -1.8460 0.8546 1.0000 2.0000 1.0000 1.0000 -1.9175 0.9319 1.0000 -2.0000 1.0000 1.0000 -1.9598 0.9604 1.0000 -2.0000 1.0000 1.0000 -1.9885 0.9889

g = 7.1024e-06

figure

freqz(b,a, 2^18, 1000)

set(subplot(2,1,1), 'XLim',[0 50])

set(subplot(2,1,2), 'XLim',[0 50])

Check = nnz(isnan(chan_data))

Check = 0

L = numel(chan_data);

L = 7000

t = linspace(0, L-1, L)/Fs;

figure

plot(t, chan_data)

grid

Fn = Fs/2;

NFFT = 2^nextpow2(L)

NFFT = 8192

FTcd = fft((chan_data-mean(chan_data)).*hann(L), NFFT)/L;

Fv = linspace(0, 1, NFFT/2+1)*Fn;

Iv = 1:numel(Fv);

figure

plot(Fv, abs(FTcd(Iv))*2)

grid

chan_filtered = filtfilt(b, a, chan_data)

chan_filtered = 7000×1

NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN

chan_filtered = filtfilt(sos, g, chan_data);

figure

plot(t, chan_filtered)

grid

chan_filtered = bandpass(chan_data, [3 20], Fs, 'ImpulseResponse','iir');

figure

plot(t, chan_filtered)

grid

The second-order-section implementation of the same Butterworth filter then works, and produces essentially the same result as an elliptic bandpass filter with the same essential parameters.

.

2 个评论
显示无隐藏无

Aboltabol 2023-5-11

编辑：Aboltabol 2023-5-11

在 MATLAB Online 中打开

Okay, so, on running this code-block:

load('Channel_Sim.mat') % Sampling frequency is 1000 Hz.
[b,a] = butter(4, [3 20] ./ (1000/2));
[z,p,k] = butter(4, [3 20] ./ (1000/2)); 
[sos,g] = zp2sos(z,p,k);
chan_filtered_1 = filtfilt(b, a, chan_data);
chan_filtered_2 = filtfilt(sos, g, chan_data);

chan_filtered_1 is NaN but chan_filtered_2 is all right. This is just weird, because they are equivalent representations?! Why does using a SOS implementations skirts around the issue?

Star Strider 2023-5-11

The second-order-section implementation of a filter is always stable (at least in my experience).

The reason has to do with the way the filters are implemented and how filtfilt handles them. The relevant discussion is in the sos documentation section for filtfilt, and a more extensive discussion is in the sosfilt documentation. (I cannot determine from reading the documentation if sosfilt does the same forward-backward filtering that filtfilt does, so filtfilt still appears to be the preferred method of doing the actual filtering.)

.

请先登录，再进行评论。