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

2 次查看（过去 30 天）
Darcy Cordell
Darcy Cordell 2023-2-10
评论： Star Strider 2023-2-14
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!

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

回答（1 个）

Star Strider
Star Strider 2023-2-11
编辑：Star Strider 2023-2-11
Ran in:
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.
.
  4 个评论
显示 2更早的评论
Darcy Cordell
Darcy Cordell 2023-2-14
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
Star Strider 2023-2-14
Ran in:
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')
.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Digital and Analog Filters 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by