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

5 views (last 30 days)
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');
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?

Answers (1)

Star Strider
Star Strider on 11 Feb 2023
Edited: Star Strider on 11 Feb 2023
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);
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.
Star Strider
Star Strider on 14 Feb 2023
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',''));
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);
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);
plot(t, data, 'DisplayName','Original')
hold on
plot(t, data_filt, 'DisplayName', 'Filtered')
hold off
NFFT = 2^nextpow2(L);
FTss = fft([data data_filt].*hann(L), NFFT);
Fv = linspace(0, 1, NFFT/2+1)*fNy;
Iv = 1:numel(Fv);
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')
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')

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!