New to filter design: Designing filters at very low frequencies (mHz) on logarithmic frequency scale

I am new to filter design, and I am trying to filter some signals to analyze the frequency content. I would like to look at low-, high-, and band-passed signals. However, I am having no end to the trouble in getting a simple butterworth filter to behave as I would like. I am sure there is some obvious misunderstanding I am having.

I want to design a filter with a passband between 10^-4 and 10^-3 Hz. Ideally, I would like the stopband corner to be very close to the passband corner in logarithmic space (e.g. 10^-4.1 and 10^-2.9).

My sample rate is 1 min (e.g. 0.0166 Hz). Here is some sample code:

fs = 1/60; %1 minute sample rate, 0.0166 Hz
fNy = fs/2; %Nyquist
Wp = 10.^[-4 -3]; %passband cutoffs
Ws = 10.^[-4.1 -2.9]; %stopband cutoffs
[n,Wn] = buttord(Wp/fNy,Ws/fNy,5,60);
[b,a] = butter(n,Wn,'bandpass');  
figure(1);
freqz(b,a,1440,fs); set(gca,'XScale','log');

This results in warnings about matrix being singular and a butterworth filter which is nonsensical since the peak of the "passband" is at -250 dB, and phase is undefined:

If I change the stopband cutoffs to be less sharp:

Ws = 10.^[-4.5 -2.5]

then I get sensible results as expected:

The filter seems very sensitive to this stopband cutoff. Now, ideally I would like the cutoff to be much sharper than what is shown in the above figure.

Can anyone explain why things are going awry with these stopband cutoff frequencies?

Does it have something to do with using logarithmic frequencies to design the filter?

Is a butterworth filter not appropriate for this task? If not, what other filter should I be using?

Thanks!

Answer 1

[Star Strider]

It is quite common for a transfer function implementation to be unstable. The solution is to use zero-pole-gain output from the filter design function, and a second-order-section implementation —

fs = 1/60; %1 minute sample rate, 0.0166 Hz
fNy = fs/2; %Nyquist
Wp = 10.^[-4 -3]; %passband cutoffs
Ws = 10.^[-4.1 -2.9]; %stopband cutoffs
[n,Wn] = buttord(Wp/fNy,Ws/fNy,5,60);
[z,p,k] = butter(n,Wn,'bandpass');
[sos,g] = zp2sos(z,p,k);

figure(1);
freqz(sos,1440,fs); set(gca,'XScale','log');

This produces a stable filter that does exactly what you requested it to do in your original code.

See the documentation on zp2sos for details. Also, always use filtfilt to do the actual filtering.

EDIT — (11 Feb 2023 at 2:02)

Added clarification.

.

[Darcy Cordell]

Thanks for the info! Can you comment on why the magnitude of the passband is ~400 dB

[Star Strider]

My pleasure!

The freqz function does not normalise its y-axis limits to the maximum of the passband, other than for a lowpass filter.

One way of getting around that is —

fs = 1/60; %1 minute sample rate, 0.0166 Hz
fNy = fs/2; %Nyquist
Wp = 10.^[-4 -3]; %passband cutoffs
Ws = 10.^[-4.1 -2.9]; %stopband cutoffs
[n,Wn] = buttord(Wp/fNy,Ws/fNy,5,60);
[z,p,k] = butter(n,Wn,'bandpass');
[sos,g] = zp2sos(z,p,k);

figure(1);
freqz(sos,1440,fs)
sp211 = subplot(2,1,1);                                                 % Handle To The Magnitude Subplot
sp212 = subplot(2,1,2);                                                 % Handle To The Phase Subplot
h = sp211.Children.YData;
sp211.Children.YData = h - max(h(isfinite(h)));                         % Correct Magnitude Subplot
sp211.XScale = 'log';
sp212.XScale = 'log';

This is a bit cumbersome, however it gives the correct plot. Only ‘sp211’ is actually required to get the magnitude normalisation correct, however you want to scale both subplot x-axes as 'log' so I retrieved both of them here to do that to both subplot axes.

EDIT — (13 Feb 2023 at 20:05)

Corrected typographical eror.

.

[Darcy Cordell]

Continuing with this, when I take the sos and g outputs from zp2sos and try to filter my signal, I get this:

fs = 1/60; %1 minute sample rate, 0.0166 Hz
fNy = fs/2; %Nyquist
Wp = 10.^[-4 -3]; %passband cutoffs
Ws = 10.^[-4.1 -2.9]; %stopband cutoffs
[n,Wn] = buttord(Wp/fNy,Ws/fNy,5,60);
[z,p,k] = butter(n,Wn,'bandpass');
[sos,g] = zp2sos(z,p,k);

%Make some data:
load('data_example.mat'); %see attachment

data_filt = filtfilt(sos,g,data);

And if I fft the filtered result, I get this:

Once again, if I change the limits on the stopband cutoffs to:

Ws = 10.^[-4.2 -2.8]; %stopband cutoffs

Then it gives sensible results (although I'm still not sure if they are "correct"):

Part of the issue is that I want to cut the signal more than is currently being cut. It seems like the stopband cutoff is not being applied correctly since at e.g. 10^-6 Hz, the "bandpass filtered" signal still has a >10^5 amplitude.

[Star Strider]

It appears that the filter is operation correctly.

What do you want it to do?

Also, if you want steeper cutoffs, a Buterworth design may not be the best option. Consider using an elliptic filter instead, and with greater stopband attenuation. (The relevant functions are ellipord and ellip, with the rest of the code unchanged otherwise.)

Try something like this —

LD = load(websave('data_example','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1295845/data_example.mat'));
data = LD.data;

fs = 1/60; %1 minute sample rate, 0.0166 Hz
fNy = fs/2; %Nyquist
Wp = 10.^[-4 -3]; %passband cutoffs
Ws = 10.^[-4.1 -2.9]; %stopband cutoffs
Rp =   1;                                    % PAssband Ripple (Attenuation)
Rs =  90;                                    % Stopband Ripple (Attenuation)
[n,Wp] = ellipord(Wp/fNy,Ws/fNy,Rp,Rs);
[z,p,k] = ellip(n,Rp,Rs,Wp,'bandpass');
[sos,g] = zp2sos(z,p,k);

figure
freqz(sos, 2^16, fs)
sp211 = subplot(2,1,1);                                                 % Handle To The Magnitude Subplot
sp212 = subplot(2,1,2);                                                 % Handle To The Phase Subplot
h = sp211.Children.YData;
sp211.Children.YData = h - max(h(isfinite(h)));                         % Correct Magnitude Subplot
sp211.XScale = 'log';
sp212.XScale = 'log';

%Make some data:
% load('data_example.mat'); %see attachment

L = numel(data);
t = linspace(0, L-1, L)/fs;

data_filt = filtfilt(sos,g,data);

figure
plot(t, data, 'DisplayName','Original')
hold on
plot(t, data_filt, 'DisplayName', 'Filtered')
hold off
legend('Location','best')
xlabel('Time')
ylabel('Amplitude')

NFFT = 2^nextpow2(L);
FTss = fft([data data_filt].*hann(L), NFFT);
Fv = linspace(0, 1, NFFT/2+1)*fNy;
Iv = 1:numel(Fv);

figure
hp{1} = semilogx(Fv, mag2db(abs(FTss(Iv,1))*2), 'DisplayName','Original');
hold on
hp{2} = plot(Fv, mag2db(abs(FTss(Iv,2))*2), 'DisplayName', 'Filtered');
hold off
% legend([hp{:}],'Location','best')
xlabel('Frequency')
ylabel('Magnitude (dB)')
title('Elliptic Filter')
xWp = xline(Wp*fNy,'-g', 'DisplayName','Passband');
xWs = xline(Ws,'-r', 'DisplayName','Stopband');
legend([hp{:} xWp(1) xWs(1)],'Location','best')

.