New to filter design: Designing filters at very low frequencies (mHz) on logarithmic frequency scale
6 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
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!
0 comentarios
Respuestas (1)
Star Strider
el 11 de Feb. de 2023
Editada: Star Strider
el 11 de Feb. de 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);
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.
EDIT — (11 Feb 2023 at 2:02)
Added clarification.
.
4 comentarios
Star Strider
el 14 de Feb. de 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','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')
.
Ver también
Categorías
Más información sobre Digital Filter Design en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!