Band pass filtering of time series data

Hello, I am trying to replicate this study that filters the time series of the H-component of the magnetic field in the period range of 10-45 seconds to see ultralow frequency variations around the time of an earthquake or volcanic event, so may I ask anyone how to use the bandpass function in MATLAB right in this case?

Here's my silly attempt

fs = 1/60; %since the data is per minute

bandpass(T.DAVH,[0.022 0.1],fs); % period range of 10-45 seconds to frequency in Hz

but it returns

Warning: Specified passband frequency is beyond the Nyquist range so signal has been
filtered with an allstop filter. 
> In highpass>designFilter (line 129)
In highpass (line 97)
In bandpass>performHighpassFiltering (line 149)
In bandpass (line 116) 

Attached the data below where T.DAVH is the H-component of the magnetic field measured in Davao station.

The description of the data is as follows:

Format IAGA-2002 |

Source of Data Kyushu University (KU) |

Station Name Davao |

IAGA CODE DAV (KU code) |

Geodetic Latitude 007.000 |

Geodetic Longitude 125.400 |

Elevation 8888.88 |

Reported HDZF |

Sensor Orientation HDZ |

Digital Sampling 1 seconds |

Data Interval Type Averaged 1-minute (00:30 - 01:29) |

Data Type Provisional

Thank you to anyone who can help!

[Mathieu NOE]

hello

there is a slight problem with the sampling rate ... if I understand correctly you have one sample per minute so fs = 1/60 = 0.0167 Hz

you cannot apply a digital band pass filter which frequencies [0.022 0.1]are above the sampling rate (Nyquist theorem)

[Mathieu NOE]

according to your post , the raw data was actually acquired with a sampling rate of 1 Hz but then averaged on a 1 minute long window

Digital Sampling 1 seconds |

Data Interval Type Averaged 1-minute (00:30 - 01:29) |

what we need here is the raw data

[Lalaine Anne Asares]

Will do! Thank you so much!

Answer 1

[Star Strider]

The problem is that the sampling frequency is 1 sample/minute, making the Nyquist frequency (the highest uniquely resolvable frequency in a sampled signal) 0.5 samples/minute. The desired stopband edge of 45 seconds (0.75 samples/minute) is higher than that.

The real problem however is that with such a low samplinng frequency (1/60 Hz), it is not possible to even record anything in that range. It might be possible to confabulate data (this requires actually creating it where it did not previously exist), however that is not likely to result in any useful information. I did that here, and the Fourier transform of your signal (upsampled) demonstrates that even in the upsampled signal, there is no energy in the range you are interested in.

My analysis —

LD = load('magneticfielddata.mat');
T = LD.T

T = 1440x8 table
       DATE         TIME      DOY    DAVH      DAVD       DAVZ     DAVF         |     
    __________    ________    ___    _____    _______    ______    _____    __________

    2023-11-29    00:00:00    333    39583    -739.78    193.28    39590    {0x0 char}
    2023-11-29    00:01:00    333    39584    -739.63    193.41    39591    {0x0 char}
    2023-11-29    00:02:00    333    39584    -739.63    193.49    39592    {0x0 char}
    2023-11-29    00:03:00    333    39584    -739.75    193.48    39591    {0x0 char}
    2023-11-29    00:04:00    333    39583    -739.85    193.38    39591    {0x0 char}
    2023-11-29    00:05:00    333    39583    -739.98     193.2    39590    {0x0 char}
    2023-11-29    00:06:00    333    39582    -740.09    192.97    39590    {0x0 char}
    2023-11-29    00:07:00    333    39582    -740.34    192.81    39590    {0x0 char}
    2023-11-29    00:08:00    333    39582    -740.15    192.61    39590    {0x0 char}
    2023-11-29    00:09:00    333    39582    -740.03    192.34    39589    {0x0 char}
    2023-11-29    00:10:00    333    39581    -740.29    192.08    39588    {0x0 char}
    2023-11-29    00:11:00    333    39580     -740.3    191.79    39588    {0x0 char}
    2023-11-29    00:12:00    333    39581    -740.36    191.59    39588    {0x0 char}
    2023-11-29    00:13:00    333    39582    -740.28    191.47    39589    {0x0 char}
    2023-11-29    00:14:00    333    39582    -739.99    191.37    39590    {0x0 char}
    2023-11-29    00:15:00    333    39582    -739.92    191.29    39589    {0x0 char}

Te = T;
DateTime = Te.DATE + Te.TIME;
DateTime.Format = 'yyyy-MM-dd HH:mm:ss';
Te = addvars(Te, DateTime, Before=1)

Te = 1440x9 table
         DateTime             DATE         TIME      DOY    DAVH      DAVD       DAVZ     DAVF         |     
    ___________________    __________    ________    ___    _____    _______    ______    _____    __________

    2023-11-29 00:00:00    2023-11-29    00:00:00    333    39583    -739.78    193.28    39590    {0x0 char}
    2023-11-29 00:01:00    2023-11-29    00:01:00    333    39584    -739.63    193.41    39591    {0x0 char}
    2023-11-29 00:02:00    2023-11-29    00:02:00    333    39584    -739.63    193.49    39592    {0x0 char}
    2023-11-29 00:03:00    2023-11-29    00:03:00    333    39584    -739.75    193.48    39591    {0x0 char}
    2023-11-29 00:04:00    2023-11-29    00:04:00    333    39583    -739.85    193.38    39591    {0x0 char}
    2023-11-29 00:05:00    2023-11-29    00:05:00    333    39583    -739.98     193.2    39590    {0x0 char}
    2023-11-29 00:06:00    2023-11-29    00:06:00    333    39582    -740.09    192.97    39590    {0x0 char}
    2023-11-29 00:07:00    2023-11-29    00:07:00    333    39582    -740.34    192.81    39590    {0x0 char}
    2023-11-29 00:08:00    2023-11-29    00:08:00    333    39582    -740.15    192.61    39590    {0x0 char}
    2023-11-29 00:09:00    2023-11-29    00:09:00    333    39582    -740.03    192.34    39589    {0x0 char}
    2023-11-29 00:10:00    2023-11-29    00:10:00    333    39581    -740.29    192.08    39588    {0x0 char}
    2023-11-29 00:11:00    2023-11-29    00:11:00    333    39580     -740.3    191.79    39588    {0x0 char}
    2023-11-29 00:12:00    2023-11-29    00:12:00    333    39581    -740.36    191.59    39588    {0x0 char}
    2023-11-29 00:13:00    2023-11-29    00:13:00    333    39582    -740.28    191.47    39589    {0x0 char}
    2023-11-29 00:14:00    2023-11-29    00:14:00    333    39582    -739.99    191.37    39590    {0x0 char}
    2023-11-29 00:15:00    2023-11-29    00:15:00    333    39582    -739.92    191.29    39589    {0x0 char}

% Te.Properties.VariableNames
Te = removevars(Te, ["DATE","TIME","|"]);
TTe = table2timetable(Te);
TTe = retime(TTe, 'secondly', 'pchip')

TTe = 86341x5 timetable
         DateTime          DOY    DAVH      DAVD       DAVZ     DAVF 
    ___________________    ___    _____    _______    ______    _____

    2023-11-29 00:00:00    333    39583    -739.78    193.28    39590
    2023-11-29 00:00:01    333    39583    -739.78    193.28    39590
    2023-11-29 00:00:02    333    39583    -739.77    193.29    39590
    2023-11-29 00:00:03    333    39583    -739.77    193.29    39590
    2023-11-29 00:00:04    333    39583    -739.77    193.29    39590
    2023-11-29 00:00:05    333    39583    -739.76    193.29    39590
    2023-11-29 00:00:06    333    39583    -739.76     193.3    39590
    2023-11-29 00:00:07    333    39583    -739.75     193.3    39590
    2023-11-29 00:00:08    333    39583    -739.75     193.3    39590
    2023-11-29 00:00:09    333    39583    -739.75     193.3    39590
    2023-11-29 00:00:10    333    39583    -739.74    193.31    39590
    2023-11-29 00:00:11    333    39583    -739.74    193.31    39590
    2023-11-29 00:00:12    333    39583    -739.74    193.31    39590
    2023-11-29 00:00:13    333    39583    -739.73    193.31    39590
    2023-11-29 00:00:14    333    39583    -739.73    193.32    39590
    2023-11-29 00:00:15    333    39583    -739.72    193.32    39590

ffttime = 0:height(TTe)-1;                                  % Time In Seconds
[FTs1,Fv] = FFT1(TTe.DAVH,ffttime(:));

figure
semilogx(Fv, abs(FTs1)*2)
grid
xlim('tight')
xline([1/45 1/10], '--r', 'Passband')
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('Fourier Transform of Upsampled DAVH')

Fs = 1;                                                 % Samples/Second (Hz)
Fn = Fs/2                                               % Nyquist Frequency

Fn = 0.5000

passband = [1/45 1/10]                                  % Passband Frequencies

passband = 1×2

0.0222 0.1000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

DAVH_filtered = bandpass(TTe.DAVH,passband,Fs); % period range of 10-45 seconds to frequency in Hz

figure
tiledlayout(2,1)
nexttile
plot(TTe.DateTime, TTe.DAVH)
grid
title('Original Signal')
nexttile
plot(TTe.DateTime, DAVH_filtered)
grid
title('Bandpass-Filtered Signal')

function [FTs1,Fv] = FFT1(s,t)
% Arguments:
%    s: Signal Vector Or Matrix
%    t: Associated Time Vector 
    t = t(:);
    L = numel(t);
    if size(s,2) == L
        s = s.';
    end
    Fs = 1/mean(diff(t));
    Fn = Fs/2;
    NFFT = 2^nextpow2(L);
    FTs = fft((s - mean(s)) .* hann(L).*ones(1,size(s,2)), NFFT)/sum(hann(L));
    Fv = Fs*(0:(NFFT/2))/NFFT;
%     Fv = linspace(0, 1, NFFT/2+1)*Fn;
    Iv = 1:numel(Fv);
    Fv = Fv(:);
    FTs1 = FTs(Iv,:);
end

There is a brief filter transient in the beginning that results from the mean of the original signal not being zero, however there is nothing after that because there is no energy in this signal in the frequency range you are interested in. The sampling freequency of the instrumentation completely failed to capture it. If you can find a recording with a sampling frequency of 1 Hz, then you could filter it to see the signal in the region of interest.

.

[Lalaine Anne Asares]

Oh, I get it now, how to bandpass filter, that is.

I was wrong to apply it in the 1-minute corrected data instead of the raw 1-second data.

Thank you so much for being detailed even though the data I attached won't bring out any useful information.

[Star Strider]

As always, my pleasure!

Use data sampled at a higher sampling frequency (at least 1 Hz, faster is always better) and you should be able to get the information you want.

[Lalaine Anne Asares]

Did it with 1 Hz sampling frequency and used 'ImpulseResponse','iir' name-value pair as per literature and looks promising so far.

The gray dashed line marks a 6.9 earthquake in the Philippines. There's still a lot to do with the data, but eyeballing the day before the earthquake, there's a large amplitude variation (max-min) which must be nature's sign that there's an upcoming earthquake.

Thank you once again!

[Star Strider]

As always, my pleasure!

If you can somehow predict earthquakes, that would be worthy of a Nobel Prize! I grew up in southern California (western U.S.), and admit to being scared of them.

I forgot to mention my preference for using ImpulseResponse='iir' since that designs a very efficient elliptic filter. (My apologies for that oversight!)

Also, you can get the amplitude ‘outline’ of the filtered waveform using the envelope function and the 'peak' option. This can be useful for estimating the beginning and ending times, and the widths at different times and amplitudes, and fitting them to nonlinear functions.

To get the amplitudes, locations, and widths (such as full-width-half-maximum) of the peaks, use the findpeaks function. Negate the filtered waveform and use findpeaks to get the ‘valleys’ and their parameters.