Band-pass filter extraction.

Dear all, How are you?

Currently, I’m dealing with a band-pass filter that extracts a specific frequency band from a synthetic time series. Unfortunately, my band-pass filter is not extracting the original series as the following code shows. I would appreciate any suggestions on how to make the filter process more effective.

clear all,close all;
Inertial_period = 20.215;
K1_period = 24;   % K1 period (hours)
M2_period = 12;
Inertial_amplitude = 1.0;    
K1_amplitude = 0.3;          
M2_amplitude = 0.5;         
Inertial_phase = 0;          % Example phase (radians)
K1_phase = pi/4;             
M2_phase = pi/2;   

% Time Vector (1 Month)
Fs = 1;                 % Sampling frequency (1 sample per hour)
duration_hours = 510; % 510 Hours
t = 0:1/Fs:duration_hours;

% Generate Components
I = Inertial_amplitude * cos(2*pi*t/Inertial_period + M2_phase);
K1 = K1_amplitude * cos(2*pi*t/K1_period + K1_phase);
M2 = M2_amplitude * cos(2*pi*t/M2_period + M2_phase);

u = K1+M2+I; %Time Serie
clearvars -except u I;

Now, I built a band-pass filter to extract the original signal "I" from the time series 'u'.

Fs = 1; % 1 sample per hour
In = 20.215;
low_cutoff = 1/(In+1);  % I + 1hour
high_cutoff = 1/(In-1); % I - 1Hour
% Normalize the cutoff frequencies by the Nyquist frequency (Fs/2)
Wn = [low_cutoff, high_cutoff] / (Fs/2);
% Design the Butterworth bandpass filter
[b, a] = butter(2, Wn, 'bandpass');
clear Fs low_cutoff_I high_cutoff_I Wn_I

fu = filtfilt(b,a,u);

Clearly, the band-pass filter with a central frequency +/- 1 hour is more accurate than the same band-pass filter with +/- 0.2 hours, but it is not identical to the original signal (orange). Is there any way I can better extract the frequency so the filtered signal will be almost identical to the original? Thanks to all; I would appreciate any suggestions.

[Star Strider]

The filter removes signal energy from the original signal. That the amplitude of the filtered signal is less than the amplitude of the original signal is to be expected. The design frequency of the filter passband appears to be about 49.5.

Fs = 1; % 1 sample per hour
In = 20.215;
low_cutoff = 1/(In+1);  % I + 1hour
high_cutoff = 1/(In-1); % I - 1Hour
% Normalize the cutoff frequencies by the Nyquist frequency (Fs/2)
Wn = [low_cutoff, high_cutoff] / (Fs/2);
% Design the Butterworth bandpass filter
[b, a] = butter(2, Wn, 'bandpass');
% clear Fs low_cutoff_I high_cutoff_I Wn_I

[h,fv] = freqz(b,a,2^16,Fs);

[Hmax,idx] = findpeaks(abs(h));
FreqMaxMag = fv(idx)/max(fv)*500

FreqMaxMag = 49.5308

figure
freqz(b,a,2^16,Fs);
xline(subplot(2,1,1),fv(idx)/max(fv)*500)
xline(subplot(2,1,2),fv(idx)/max(fv)*500)
set(subplot(2,1,1), 'XGrid','on', 'YGrid','on')
set(subplot(2,1,2), 'XGrid','on', 'YGrid','on')

.

Answer 1

[William Rose]

@Claudio Iturra,

Use a lower low frequency and a higher high frequency when specifying the bandpass filter, and consider using a higher order filter.

Let's consider a simple lowpass filter with filtfilt(). If the filter designer specifies a Butterworth with a cutoff frequency f0 (hoping for 3 dB attenuation at f0), the Butterworth will be 3 dB down at f=f0. But with filtfilt(), the attenuation will be 6 dB at f0, since the signal gets passed through the filter twice. A solution is to specify the cutoff of the initial lowpass Butterworth as 1.25*f0. For a second order Butterworth lowpass filter, this results in 1.5 dB attenuation at f0. After the signal is filtered twice with filtfilt(), there will be 3 dB attenuation at f0, which is what the designer wants.

Since you are using a bandpass, the details will differ somewhat, but you get the idea: "widen" the bandpass, and use a higher order filter. This will result in the passband amplitude being closer to unity. This will extract the "I" (inertial) signal with an amplitude closer to its orignal amplitude.

[William Rose]

@Claudio Iturra,

My comment above is based on the signal where the filter cutoff perfiods are +- 0.2 hours. IN that case, the amplitude of the passband is less than unity, and my suiggesitons will bring it closer to unity.

When you used filtercutoff periods +-1 hour (wider pass band, as I recommended), the passband amplitude IS close to unity, which is why the extracted signal is very close to equal to the desired cmponent, in the middle of the siognal. However, it is not so close at the edges. If you make the cutoff periods +-2 hours, the result is better. See figure below. Filter 1 uses +- 1 hour, fiter 2 uses +- 2 hours. Both filters are second order, before using filtfilt().

You will never get the initial and final transients perfect, especially with bandpass filters. In fact, if you use too high an order, and too narrow a band, you can get significant ringing at the start and end (when using filtfilt()), that takes a while to die out. This is an unavoidable consequence of using a sharp and narrow bandpass filter.

Answer 2

[Paul]

Hi Claudio,

Start with the original code

Inertial_period = 20.215;
K1_period = 24;   % K1 period (hours)
M2_period = 12;
Inertial_amplitude = 1.0;    
K1_amplitude = 0.3;          
M2_amplitude = 0.5;         
Inertial_phase = 0;          % Example phase (radians)
K1_phase = pi/4;             
M2_phase = pi/2;   
% Time Vector (1 Month)
Fs = 1;                 % Sampling frequency (1 sample per hour)
duration_hours = 510; % 510 Hours
t = 0:1/Fs:duration_hours;
% Generate Components
I = Inertial_amplitude * cos(2*pi*t/Inertial_period + M2_phase);
K1 = K1_amplitude * cos(2*pi*t/K1_period + K1_phase);
M2 = M2_amplitude * cos(2*pi*t/M2_period + M2_phase);
u = K1+M2+I; %Time Serie
%clearvars -except u I;

%Now, I built a band-pass filter to extract the original signal "I" from the time series 'u'.
Fs = 1; % 1 sample per hour
In = 20.215;
low_cutoff = 1/(In+1);  % I + 1hour
high_cutoff = 1/(In-1); % I - 1Hour
% Normalize the cutoff frequencies by the Nyquist frequency (Fs/2)
Wn = [low_cutoff, high_cutoff] / (Fs/2);
% Design the Butterworth bandpass filter
[b, a] = butter(2, Wn, 'bandpass');
%clear Fs low_cutoff_I high_cutoff_I Wn_I

At this point, it's probably easier to see what's going on if we apply the fitler to I instead of u. Let's filter forward only

fu = filter(b,a,I);
figure
plot(t,I,t,fu)

The bandpass effect of the filter is a steady-state concept. But, as we can see from the plot, the filter has a transient response of roughly 200 seconds before the output reaches steady state.

Now, when we use filtfilt to do forward/backard filtering, we're kind of, sort of, essentially, approximately also applying the filter to the red curve above, but starting from the right and going to the left. So we'd expect that same 200 second transient, but going from right to left,

fu = filtfilt(b,a,I);
figure
plot(t,I,t,fu)

which is basically what happens. There's a very small window where the response has reached steady state going in both directions.

So I believe the effect you're observing is due to the slow transient reponse of the filter relative to the time span of the data. Let's try with more data.

t = 0:1/Fs:3*duration_hours;
% Generate Components
I = Inertial_amplitude * cos(2*pi*t/Inertial_period + M2_phase);
figure
plot(t,I,t,filtfilt(b,a,I))

The effect of the transient response is still there, of course, but now the red curve looks more like what's probably expected.

If more data can't be used, then I guess one could try to design a bandpass filter with a much faster transient response.