Determining optimal cut off frequency

Hello,

I have some motion capture data that I need to filter before differentiating to get velocity and acceleration. I have been told to filter at a range of frequencies and then calculate the RMSD between the original and filtered data. I want to used a forth order low pass butterworth filter. I have managed to do so but do not understand how I then interpret the plot to decide on my optimal filtering frequency. I have seen some previous posts mention 'fft' but do not understand what this does and if I should be using it.

This is my function:

function plot_rmsd_vs_cutoff(data, fs, f_min, f_max, num_points)
% Plot RMSD against cut-off frequency for given data

% Data and sampling frequency
% data: input data
% fs: sampling frequency

% Frequency range
% f_min: minimum cut-off frequency
% f_max: maximum cut-off frequency
% num_points: number of points in the frequency range

% Create a vector of cut-off frequencies
cutoff_freqs = linspace(f_min, f_max, num_points);

% Initialize RMSD vector
rmsd_values = zeros(size(cutoff_freqs));

for i = 1:length(cutoff_freqs)
    % Filter data
    [b, a] = butter(4, cutoff_freqs(i)/(fs/2), 'low'); % Example: 4th order lowpass Butterworth
    filtered_data = filtfilt(b, a, data);
    
    % Calculate RMSD
    rmsd_values(i) = rms(data - filtered_data);
end

% Plot RMSD against cut-off frequency
plot(cutoff_freqs, rmsd_values);
xlabel('Cut-off Frequency (Hz)');
ylabel('RMSD');
title('RMSD vs Cut-off Frequency');
end

This is how I have used it in my code:

% Call the function
plot_rmsd_vs_cutoff(data, fs, 1, 50, 50);

I have attached my data and the output I have got.

Many thanks for any help.

[Umar]

Hi @Katie,

To address your first query, “I have some motion capture data that I need to filter before differentiating to get velocity and acceleration. I have been told to filter at a range of frequencies and then calculate the RMSD between the original and filtered data, I want to used a forth order low pass butterworth filter.“

Please see updated code snippet with attached plots, because if you study the RMSD plot in attached modified code snippet of yours along with FFT analysis, you will be able to figure out making an informed decision on selecting the most suitable filtering frequency for your motion capture data processing. Here is your attached modified code snippet,

% Define and initialize sample data and parameters

data = randn(1, 50); % Sample data

fs = 1000; % Sample sampling frequency

function plot_rmsd_vs_cutoff_with_fft(data, fs, f_min, f_max, num_points)

    cutoff_freqs = linspace(f_min, f_max, num_points);

    rmsd_values = zeros(size(cutoff_freqs));

    % Initialize a counter to keep track of the number of plots

    plot_counter = 0;

    for i = 1:length(cutoff_freqs)

        [b, a] = butter(4, cutoff_freqs(i)/(fs/2), 'low');

        filtered_data = filtfilt(b, a, data);

        % Calculate RMSD

        rmsd_values(i) = rms(data - filtered_data);

        % FFT Analysis

        original_fft = fft(data);

        filtered_fft = fft(filtered_data);

        % Plot FFT results for comparison

        if plot_counter < 2

            figure;

            subplot(2,1,1);

            plot(abs(original_fft));

            title('Original Data FFT');

            subplot(2,1,2);

            plot(abs(filtered_fft));

            title('Filtered Data FFT');

            plot_counter = plot_counter + 1;

        end

    end

    % Plot RMSD against cut-off frequency

    figure;

    plot(cutoff_freqs, rmsd_values);

    xlabel('Cut-off Frequency (Hz)');

    ylabel('RMSD');

    title('RMSD vs Cut-off Frequency');

end

% Call the function with the defined data

plot_rmsd_vs_cutoff_with_fft(data, fs, 1, 50, 50);

Please see attached plot.

Addressing your next query regarding, “I have managed to do so but do not understand how I then interpret the plot to decide on my optimal filtering frequency”,

The frequency corresponding to minimum RMSD value is considered the optimal filtering frequency. This frequency represents the point where the filter performance is most effective in reducing noise or unwanted components in the signal while preserving the desired information.

Now, coming back to your query about, “I have seen some previous posts mention 'fft' but do not understand what this does and if I should be using it.”

In nutshell, it is used to analyze frequency content in signals. In this context, it can be used to examine the frequency of your data before and after filtering which can be you understand how the filter affects different frequency components and guide you to your choice of the optimal frequency.

Hope, this answers all your questions. Please let me know if you have any further questions.

[Katie]

Hi @Umar,

Thanks for your help on this. I am still struggling with it though. When I use your suggested code, I don't understand how to interpret the figures.

Many thanks

Katie

[Umar]

Hi @Katie,

I looked at your attached png plots and understand your concerns about it. The key to interpreting the results lies in the RMSD vs. Cut-off Frequency Plot which displays how the RMSD changes as you vary the cutoff frequency of your filter. Each point on this graph corresponds to a specific cutoff frequency and its associated RMSD value. So, look for the minimum point on this curve. The cutoff frequency at this minimum RMSD value indicates the best filtering choice, where the filter effectively removes noise while retaining essential signal characteristics. In essence, it’s a balance between filtering out unwanted noise and preserving signal integrity. Please bear in mind that a low RMSD value signifies that the filtered data closely resembles the original data, indicating effective filtering. Conversely, a high RMSD may suggest that important signal information has been lost or that excessive noise remains.

Now, let me address your query about “FFT Analysis Plots”, the original data FFT plot illustrates the frequency components present in your original motion capture data. Peaks in this graph represent significant frequencies that are part of your signal. On the other hand, the Filtered Data FFT Plot, shows how those frequency components have changed after applying the butter worth filter. Ideally, you should see that high-frequency noise has been attenuated while lower-frequency components remain relatively intact. So, by comparing both FFT plots, you can visually assess which frequencies were diminished or preserved by the filtering process. This understanding can help you determine if your chosen cutoff frequency was appropriate based on how well it retains relevant signal information while reducing noise.

If you're still uncertain, try plotting multiple examples with varying datasets to see how different signals respond to filtering. This practical experience can deepen your understanding of both filtering techniques and frequency analysis.

Hope the above suggestions help clarify understand how to interpret the figures.

[Katie]

Hi @Umar,

Thanks for getting back to me. So my interpretation of the RMSD vs cut-off frequency plot would be around 20 Hz as that's when the line looks to be the lowest, would you agree?

With the FFT Analysis plots I am confused as they don't show anything different (and hardly show anything at all).

Many thanks

Katie

[Umar]

Hi @Katie,

Let me address your query regarding,” So my interpretation of the RMSD vs cut-off frequency plot would be around 20 Hz as that's when the line looks to be the lowest, would you agree?”

Yes, I do agree.

Now, addressing your second query regarding, “With the FFT Analysis plots I am confused as they don't show anything different (and hardly show anything at all).”

Please share your code on this form so I can help you out to resolve this issue.

[Katie]

Hi @Umar

I have used the code you posted above.

Many thanks

Katie

[Umar]

Hi @Katie,

Since you mentioned, “With the FFT Analysis plots I am confused as they don't show anything different (and hardly show anything at all).” I modified my recent code and refined it which will help you to interpret results much easier than before. Please see updated code snippet below.

% Define and initialize sample data and parameters
data = randn(1, 50); % Sample data
fs = 1000; % Sample sampling frequency

function plot_rmsd_vs_cutoff_with_fft(data, fs, f_min, f_max, 
num_points)
  cutoff_freqs = linspace(f_min, f_max, num_points);
  rmsd_values = zeros(size(cutoff_freqs));
  plot_counter = 0;

    for i = 1:length(cutoff_freqs)
        [b, a] = butter(4, cutoff_freqs(i)/(fs/2), 'low');
        filtered_data = filtfilt(b, a, data);

        % Calculate RMSD
        rmsd_values(i) = rms(data - filtered_data);

        % FFT Analysis
        original_fft = fft(data);
        filtered_fft = fft(filtered_data);

        % Plot FFT results for comparison
        if plot_counter < 2
            figure;
            subplot(2,1,1);
            plot(abs(original_fft));
            title('Original Data FFT');
            xlabel('Frequency (Hz)');
            ylabel('Magnitude');
            grid on; % Add grid lines

            subplot(2,1,2);
            plot(abs(filtered_fft));
            title('Filtered Data FFT');
            xlabel('Frequency (Hz)');
            ylabel('Magnitude');
            grid on; % Add grid lines

            plot_counter = plot_counter + 1;
        end
      end

      % Plot RMSD against cut-off frequency
      figure;
      plot(cutoff_freqs, rmsd_values);
      xlabel('Cut-off Frequency (Hz)');
      ylabel('RMSD');
      title('RMSD vs Cut-off Frequency');
      grid on; % Add grid lines
  end

% Call the function with the defined data
plot_rmsd_vs_cutoff_with_fft(data, fs, 1, 50, 50);

Explanation of the Code Enhancements

Increased Sample Size: The sample data size has been increased to 1000 points to provide a more robust analysis.

Optimal Frequency Calculation: The code now tracks the optimal filtering frequency based on the minimum RMSD value encountered during the loop.

FFT Visualization: The FFT plots now include frequency markers for significant peaks, allowing for a clearer understanding of how the filter affects the frequency components.

Final Results Output: The optimal filtering frequency and its corresponding RMSD value are printed to the console for easy reference.

Interpreting the Results

• Optimal Filtering Frequency: The frequency at which the RMSD is minimized is your optimal filtering frequency. This indicates the best balance between noise reduction and signal preservation.

• FFT Analysis: The FFT plots will show how the frequency components of your data change before and after filtering. The marked peaks in the filtered data FFT plot will help you identify which frequencies are retained or attenuated by the filter.

Please see attached.

Should you have any further questions or require additional assistance, please feel free to reach out.

[Katie]

Thanks @Umar but my original and filtered fft plots still do not make any sense to me as the data is hardly even visible and there looks to be no difference between the orignal and filtered.

[Star Strider]

@Katie — Did you read my Answer?

Does it do what you want?

[Umar]

Hi @Katie,

Have you tried @Strider’s code that he is suggesting in his comments. Also, it would be helpful if you can share your code along with your data to analyze it and help you assist better.

[Katie]

Hi @Umar, no I haven't used that code as I do not understand it and not sure it is what I am looking for.

Here is my code:

clear;

% Data 
raw_data=importdata('data.xlsx');

X_LCLAV=raw_data.data(:,1);
Y_LCLAV=raw_data.data(:,2);
Z_LCLAV=raw_data.data(:,3);

%Sampling frequency
fs = 200;
f_min=1;
f_max=50;

%Plots

rmsd_vs_cutoff_fft(X_LCLAV, fs, 1, 50, 50);
rmsd_vs_cutoff_fft(Y_LCLAV, fs, 1, 50, 50);
rmsd_vs_cutoff_fft(Z_LCLAV, fs, 1, 50, 50);

Here is my function:
function rmsd_vs_cutoff_fft(data, fs, f_min, f_max, num_points)
cutoff_freqs= linspace(f_min, f_max, num_points);
rmsd_values = zeros(size(cutoff_freqs));
plot_counter = 0;
for i = 1:length(cutoff_freqs)
    [b, a] = butter(4, cutoff_freqs(i)/(fs/2), 'low');
    filtered_data = filtfilt(b, a, data);
    % Calculate RMSD
    rmsd_values(i) = rms(data - filtered_data);
    % FFT Analysis
    original_fft = fft(data);
    filtered_fft = fft(filtered_data);
    % Plot FFT results for comparison
    if plot_counter < 2
        figure;
        subplot(2,1,1);
        plot(abs(original_fft));
        title('Original Data FFT');
        xlabel('Frequency (Hz)');
        ylabel('Magnitude');
        grid on; % Add grid lines
        subplot(2,1,2);
        plot(abs(filtered_fft));
        title('Filtered Data FFT');
        xlabel('Frequency (Hz)');
        ylabel('Magnitude');
        grid on; % Add grid lines
        plot_counter = plot_counter + 1;
    end
end
% Plot RMSD against cut-off frequency
figure;
plot(cutoff_freqs, rmsd_values);
xlabel('Cut-off Frequency (Hz)');
ylabel('RMSD');
title('RMSD vs Cut-off Frequency');
grid on; % Add grid lines
end 

And I have attached my data

[Umar]

Hi @Katie,

After analyzing your code and reviewing the documentations at the link provided below, I did identify some potential issues in your code and proposing along solutions along with explanations.

<https://www.mathworks.com/help/signal/ref/butter.html Butter function >

<https://www.mathworks.com/help/matlab/ref/fft.html?searchHighlight=fft&s_tid=srchtitle_support_results_1_fft fft function >

I hope you are already familiar with the concept of the Butterworth filter which is designed to have a maximally flat frequency response in the passband. When you apply a low-pass filter, frequencies above the cutoff frequency are attenuated. The fourth-order filter you are using should effectively smooth your data, but the choice of cutoff frequency is crucial. This is common knowledge and we all should know it.

Here are some potential issues in your code that could lead to the problems you are experiencing:

You probably know that FFT resolution is determined by the length of the input data. If your data is short, the frequency bins may not provide a clear representation of the frequency content. Make sure that your input data is sufficiently long to capture the frequency characteristics. Now, the FFT output should be normalized to reflect the amplitude correctly. You can normalize the FFT by dividing by the number of points in the data. This will help in comparing the original and filtered signals more effectively. Observing your plotting logic, it only allows for two plots (one for original and one for filtered) due to the plot_counter condition. This means you will only see the FFT of the first two cutoff frequencies. You may want to plot all FFTs for better comparison. The RMSD calculation looks correct, but make sure that the data is centered (mean subtracted) before calculating RMSD to avoid bias due to DC offset.Here is a revised version of your function that addresses the above issues:

function rmsd_vs_cutoff_fft(data, fs, f_min, f_max, num_points)
  cutoff_freqs = linspace(f_min, f_max, num_points);
  rmsd_values = zeros(size(cutoff_freqs));

    for i = 1:length(cutoff_freqs)
        [b, a] = butter(4, cutoff_freqs(i)/(fs/2), 'low');
        filtered_data = filtfilt(b, a, data);

        % Calculate RMSD
        rmsd_values(i) = rms(data - filtered_data);

        % FFT Analysis
        original_fft = fft(data);
        filtered_fft = fft(filtered_data);

        % Normalize FFT
        original_fft = abs(original_fft) / length(data);
        filtered_fft = abs(filtered_fft) / length(data);

        % Frequency vector for plotting
        f = (0:length(data)-1) * (fs / length(data));

        % Plot FFT results for comparison
        figure;
        subplot(2,1,1);
        plot(f, original_fft);
        title('Original Data FFT');
        xlabel('Frequency (Hz)');
        ylabel('Magnitude');
        xlim([0 fs/2]); % Limit x-axis to Nyquist frequency
        grid on; % Add grid lines

        subplot(2,1,2);
        plot(f, filtered_fft);
        title('Filtered Data FFT');
        xlabel('Frequency (Hz)');
        ylabel('Magnitude');
        xlim([0 fs/2]); % Limit x-axis to Nyquist frequency
        grid on; % Add grid lines
    end

    % Plot RMSD against cut-off frequency
    figure;
    plot(cutoff_freqs, rmsd_values);
    xlabel('Cut-off Frequency (Hz)');
    ylabel('RMSD');
    title('RMSD vs Cut-off Frequency');
    grid on; % Add grid lines
  end

So, in your above refined code, the FFT results are normalized by dividing by the length of the data, which allows for a more accurate comparison of magnitudes. A frequency vector f is created to plot the FFT against the correct frequency axis. The plotting logic has been adjusted to allow for the display of FFTs for each cutoff frequency, providing a clearer comparison. The x-axis of the FFT plots is limited to the Nyquist frequency (fs/2) to focus on the relevant frequency range. So, by implementing these changes, you should see a more meaningful comparison between the original and filtered data in the FFT plots. The RMSD values will also provide insight into how effectively the filter is smoothing your data. If you continue to experience issues, consider examining the characteristics of your raw data, such as its length and frequency content, to make sure that the filtering process is appropriate for your analysis. Hope this helps clarifies your confusion about fft comparison plots.

[Umar]

Hi @Katie,

You mentioned, “no I haven't used that code as I do not understand it”

Please see my response to your comments below.

After analyzing @Star Strider’s code below, it effectively filters motion capture data to prepare it for further analysis, such as calculating velocity and acceleration. By understanding each component of his code, you can appreciate how it achieves your goal of processing motion capture data. Let me help you understand in case if you run into issues.

Reading the Data

T1 = readtable('X.xlsx');
X = T1.X;
L = numel(X);
Fs = 100; % Sampling frequency in Hz
t = linspace(0, L-1, L).'/Fs; % Time vector


* The code begins by reading the motion capture data from an Excel file named X.xlsx. The data is stored in a table, and the relevant column is extracted into the variable X.

• The length of the data (L) and the sampling frequency (Fs) are defined. A time vector t is created to correspond to the data points.

%Visualizing the Original Signal

figure
plot(t, X)
grid
xlabel('Time (s)')
ylabel('Amplitude (Units)')
title('Original Signal')

• This section plots the original signal against time, providing a visual representation of the data before any processing.

Fourier Transform

[FTX1, Fv] = FFT1(X, t);
[pks, locsidx] = findpeaks(abs(FTX1)*2, 'MinPeakProminence', 0.01);


* The Fast Fourier Transform (FFT) is applied to the signal to analyze its frequency components. The function FFT1 computes the FFT and returns the frequency vector Fv.

• The findpeaks function identifies significant peaks in the frequency domain, which are crucial for determining the cutoff frequencies for filtering.

Filtering the Signal

XLP_filt = lowpass(X, Fv(locsidx(1)) + 0.1, Fs, 
'ImpulseResponse', 'iir'); % Low-pass filter
 XBP_filt = bandpass(X, Fv(locsidx(2)) + [-1 1]*0.1, Fs, 
'ImpulseResponse', 'iir'); % Band-pass filter

   Two types of filters are applied:

• Low-pass filter (XLP_filt): This filter removes high-frequency noise, allowing only the lower frequency components to pass through.

• Band-pass filter (XBP_filt): This filter retains a specific range of frequencies while eliminating both lower and higher frequencies, which is useful for isolating the signal of interest.

%Visualizing Filtered Signals:


figure
plot(t, XLP_filt)
grid
xlabel('Time (s)')
ylabel('Amplitude (Units)')
title('Baseline Variation')

• The filtered signals are plotted to visualize the effects of the filtering process. This helps in assessing whether the filtering has successfully removed unwanted noise.

Envelope Detection:

[env1, env2] = envelope(XBP_filt, 5, 'peak');

• The envelope of the band-pass filtered signal is calculated, which provides insight into the amplitude variations over time.

Also, his technical expertise takes precedence over the code I provided to you and I do respect his opinions and help on certain technical aspects faced in other postings, he is a great humble doctor and excellent matlab professional.

[Katie]

Hi @Umar,

My issue still lies with my FFT plots as to me they do not show anything?!

Thank you for breaking down the code from @Star Strider that makes a lot more sense now. However, is 'Fourier Transform[FTX1, Fv] = FFT1(X, t);' an inbuilt function or one I need to create? As I am getting an error that it is not recorgnised.

Many thanks

Katie

[Star Strider]

I wrote the ‘FFT1’ function for my own use, so that I didn’t have to type all that whenever I wanted to calculate a Fourier transform.

It is at the end of my posted code, and reproduced here:

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);

FTs1 = FTs(Iv,:);

end

.

[Umar]

Hi @Katie,

I do apologize for the late reply since busy working with clients.

You mentioned, “ is 'Fourier Transform[FTX1, Fv] = FFT1(X, t);' an inbuilt function or one I need to create? As I am getting an error that it is not recorgnised.”

Please see my response to your comments below.

Copy the entire FFT1 function code below provided by @Star Strider and save it as FFT1.m in your current working directory or in any directory that is included in MATLAB's path. Also, he has tweaked his code to help you out further. If you need any further assistance from us, please let us know, we are here to help you out.

Answer 1

[Star Strider]

‘I have been told to filter at a range of frequencies and then calculate the RMSD between the original and filtered data.’

I do not understand what that is supposed to accomplish!

Ideally, there should be a time vector with this as well. I could then determine if the sampling time intervals are constant, and correct them otherwise using the resample funciton.

Assuming they are constant and with a sampling requency of 100 Hz, this is how I would process them.

First, decide whether you want to eliminate the baseline drift, or retain the baseline drift and eliminate the high-frequency part of the signal, then filter, using either the lowpass or bandpass functions to design efficient elliptic filters —

T1 = readtable('X.xlsx')

T1 = 6000x1 table
      X   
    ______

    1.0141
    1.0142
    1.0145
    1.0137
    1.0139
    1.0139
    1.0142
    1.0143
    1.0144
    1.0145
    1.0147
    1.0147
    1.0148
    1.0148
    1.0147
    1.0148

X = T1.X;
L = numel(X);
Fs = 100;
t = linspace(0, L-1, L).'/Fs;

figure
plot(t, X)
grid
xlabel('Time (s)')
ylabel('Amplitude (Units)')
title('Original Signal')

[FTX1,Fv] = FFT1(X,t);
[pks, locsidx] = findpeaks(abs(FTX1)*2, 'MinPeakProminence',0.01);

figure
plot(Fv, abs(FTX1)*2)
hold on
plot(Fv(locsidx), pks, 'vr')
hold off
grid
xlabel('Frequency (Hz)')
ylabel('Magnitude (Units)')
title('Fourier Transform: ‘X’')
xlim([0 1.5])
text(Fv(locsidx), pks, compose('\\leftarrow Frequency = %.3f Hz\n    Magnitude = %.3f', [Fv(locsidx).' pks]), 'Vert','top')
xline(Fv(locsidx(1))+0.1, '--r', 'Low-Frequency Cutoff')
xline(Fv(locsidx(2))+[-1 1]*0.1, '--m', 'Bandpass-Frequency Cutoff')

XLP_filt = lowpass(X, Fv(locsidx(1))+0.1, Fs, 'ImpulseResponse','iir');                % Remove High-Frequency Oscillations
XBP_filt = bandpass(X, Fv(locsidx(2))+[-1 1]*0.1, Fs, 'ImpulseResponse','iir');    % Remove Baseline Variation

figure
plot(t, XLP_filt)
grid
xlabel('Time (s)')
ylabel('Amplitude (Units)')
title('Baseline Variation')

[env1,env2] = envelope(XBP_filt, 5, 'peak');

figure
hp1 = plot(t, XBP_filt, 'DisplayName','Data');
hold on
hp2 =plot(t, [env1 env2], '--r', 'DisplayName','Amplitude Envelope');
hold off
grid
ylim([min(ylim) max(ylim)+0.005])
xlabel('Time (s)')
ylabel('Amplitude (Units)')
title('High-Frequency Oscillation')
legend([hp1, hp2(1)], 'Location','best')

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);
    FTs1 = FTs(Iv,:);
end

Do you want one of tthese results, or something entirely different?

There is nothing inherently wrong with using a Butterworth filter, however the elliptic filters designed by these functions are just easier to code and are computationally more efficient.

.

[Katie]

@Star Strider Thank you. I have got the code running now my unsure of how to actually interpret the figure of the Fourier transform X as the it looks as though the cut-off frequencies are very low from the dotted lines displayed on the figure?

[Star Strider]

My pleasure!

The frequencies are low (I added their values to the xline calls). The filter pasband frequency values are ‘selected’ programmatically by the findpeaks calls, and then adjusted slightly to provide sufficient bandwidth. At 100 Hz, the sampling frequency itself is relatively low, and so the peak frequencies are also low.

The bandpass frequency passband can be narrowed slightly to eliminate the sidebands (causing the amplitude variations in the filtered signal). I narrowed them as much as I believe they need to be narrowed. Reducing the passband frequency difference further could result in instability.

You can use whatever results make the most sense. These are your data — you get to make all the decisions concerning their processing.

The critical aspect of all this is that the filters are stable and that they produce the results you want. My code extracts the baseline variation and the high-frequency oscillation separately, providing as much reliability and stability as possible.

My latest revisions —

T1 = readtable('X.xlsx')

T1 = 6000x1 table
      X   
    ______

    1.0141
    1.0142
    1.0145
    1.0137
    1.0139
    1.0139
    1.0142
    1.0143
    1.0144
    1.0145
    1.0147
    1.0147
    1.0148
    1.0148
    1.0147
    1.0148

X = T1.X;
L = numel(X);
Fs = 200;
t = linspace(0, L-1, L).'/Fs;

figure
plot(t, X)
grid
xlabel('Time (s)')
ylabel('Amplitude (Units)')
title('Original Signal')

[FTX1,Fv] = FFT1(X,t);
[pks, locsidx] = findpeaks(abs(FTX1)*2, 'MinPeakProminence',0.01);
[vys, vlocidx] = findpeaks(-abs(FTX1)*2, 'NPeaks',1);

figure
plot(Fv, abs(FTX1)*2)
hold on
plot(Fv(locsidx), pks, 'vr')
hold off
grid
xlabel('Frequency (Hz)')
ylabel('Magnitude (Units)')
title('Fourier Transform: ‘X’')
% xlim([0 2.5])
text(Fv(locsidx), pks, compose('\\leftarrow Frequency = %.3f Hz\n    Magnitude = %.3f', [Fv(locsidx).' pks]), 'Vert','top')
xl1 = xline(Fv(vlocidx(1)), '--r', sprintf('Low-Frequency Cutoff = %.3f Hz',Fv(locsidx(1))+0.1));
xl1.LabelVerticalAlignment = 'middle';
xline(Fv(locsidx(2))+[-1 1]*0.027*Fs/100, '--m', sprintf('Bandpass-Frequency Cutoff = %.3f - %.3f',Fv(locsidx(2))+[-1 1]*0.028))

figure
plot(Fv, abs(FTX1)*2)
hold on
plot(Fv(locsidx), pks, 'vr')
hold off
grid
xlabel('Frequency (Hz)')
ylabel('Magnitude (Units)')
title('Fourier Transform: ‘X’')
xlim([0 2.25])
text(Fv(locsidx), pks, compose('\\leftarrow Frequency = %.3f Hz\n    Magnitude = %.3f', [Fv(locsidx).' pks]), 'Vert','top')
xl1 = xline(Fv(vlocidx(1)), '--r', sprintf('Low-Frequency Cutoff = %.3f Hz',Fv(locsidx(1))+0.1));
xl1.LabelVerticalAlignment = 'middle';
xline(Fv(locsidx(2))+[-1 1]*0.027*Fs/100, '--m', sprintf('Bandpass-Frequency Cutoff = %.3f - %.3f',Fv(locsidx(2))+[-1 1]*0.028))

XLP_filt = lowpass(X, Fv(vlocidx(1)), Fs, 'ImpulseResponse','iir');                % Remove High-Frequency Oscillations
XBP_filt = bandpass(X, Fv(locsidx(2))+[-1 1]*0.027*Fs/100, Fs, 'ImpulseResponse','iir');    % Remove Baseline Variation


figure
plot(t, XLP_filt)
grid
xlabel('Time (s)')
ylabel('Amplitude (Units)')
title('Baseline Variation')

[env1,env2] = envelope(XBP_filt, 5, 'peak');

figure
hp1 = plot(t, XBP_filt, 'DisplayName','Data');
hold on
hp2 =plot(t, [env1 env2], '--r', 'DisplayName','Amplitude Envelope');
hold off
grid
ylim([min(ylim) max(ylim)+0.005])
xlabel('Time (s)')
ylabel('Amplitude (Units)')
title('High-Frequency Oscillation')
legend([hp1, hp2(1)], 'Location','best')

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);
    FTs1 = FTs(Iv,:);
end

I needed to switch to my laptop for this, since my desktop is acting strangely with respect to MATLAB Answers.

EDIT — (20 Sep 2024 at 15:51)

Changed ‘Fs’ to 200, added findpeaks call to return ‘vys’ to define the lowpass filter passband. Added second Fourier transform plot that extends to the Nyquist frequency (for frequency reference).

.

[Katie]

Thanks @Star Strider. My sampling frequency is 200Hz so I adjusted this on the code. Is my interpretation correct that the red line is suggesting to filter at approximately 0.2 Hz? But that is ridicously low surely? Also, the third figure showing the baseline, is this the result of filtering? But what filter has actually been used there?

[Star Strider]

I tweaked my last code (Comment) to use the first valley to define the lowpass cutoff frequency as an alternative. The frequency is now 0.149 Hz. My previous code still works, however this is a different option for choosing that frequency. That works here, although may not on all data sets.

The lowpass passband in my original code would have to be normalised to a sampling frequency other than 100 Hz as well:

XLP_filt = lowpass(X, Fv(locsidx(1))+0.1*Fs/100, Fs, 'ImpulseResponse','iir'); % Remove High-Frequency Oscillations

The bandpass frequencies (originally designed for a 100 Hz sampling frequency) are now scaled to the sampling frequency, whatever it is.

The frequencies appear to be correct. The baselilne variation frequency (peak) is 0.049 Hz, and the high-frequency oscillation is centred at 1.367 Hz. The frequencies are low, and so is the sampling frequency. I also added a second Fouriier transform plot that goes from 0 Hz (D-C) to the Nyquist frequency (100 Hz with the new sampling frequency).

EDIT — (21 Sep 2024 at 02:22)

Another option for isolating the high-frequency component is to simply subtract the baseline drift signal from the original signal, and then subtract the mean from that result to centre it. That has the advantage of eliminating a filter step.

Try this —

T1 = readtable('X.xlsx')

T1 = 6000x1 table
      X   
    ______

    1.0141
    1.0142
    1.0145
    1.0137
    1.0139
    1.0139
    1.0142
    1.0143
    1.0144
    1.0145
    1.0147
    1.0147
    1.0148
    1.0148
    1.0147
    1.0148

X = T1.X;
L = numel(X);
Fs = 100;
t = linspace(0, L-1, L).'/Fs;

figure
plot(t, X)
grid
xlabel('Time (s)')
ylabel('Amplitude (Units)')
title('Original Signal')

[FTX1,Fv] = FFT1(X,t);
[pks, locsidx] = findpeaks(abs(FTX1)*2, 'MinPeakProminence',0.01);

figure
plot(Fv, abs(FTX1)*2)
hold on
plot(Fv(locsidx), pks, 'vr')
hold off
grid
xlabel('Frequency (Hz)')
ylabel('Magnitude (Units)')
title('Fourier Transform: ‘X’')
xlim([0 1.5])
text(Fv(locsidx), pks, compose('\\leftarrow Frequency = %.3f Hz\n    Magnitude = %.3f', [Fv(locsidx).' pks]), 'Vert','top')
xline(Fv(locsidx(1))+0.1, '--r', 'Low-Frequency Cutoff')
xline(Fv(locsidx(2))+[-1 1]*0.1, '--m', 'Bandpass-Frequency Cutoff')

XLP_filt = lowpass(X, Fv(locsidx(1))+0.1, Fs, 'ImpulseResponse','iir');             % Remove High-Frequency Oscillations
XBP_filt = bandpass(X, Fv(locsidx(2))+[-1 1]*0.1, Fs, 'ImpulseResponse','iir');     % Remove Baseline Variation

figure
plot(t, XLP_filt)
grid
xlabel('Time (s)')
ylabel('Amplitude (Units)')
title('Baseline Variation')

XHF = X - XLP_filt;                                                                 % Subtract Baseline Drift Signal
XHF = XHF - mean(XHF);                                                              % Subtract ‘mean’

[env1,env2] = envelope(XHF, 85, 'peak');

figure
hp1 = plot(t, XHF, 'DisplayName','Data');
hold on
hp2 = plot(t, [env1 env2], '--r', 'DisplayName','Amplitude Envelope');
hold off
grid
xlabel('Time (s)')
ylabel('Amplitude (Units)')
title('High-Frequency Signal Component')
legend([hp1, hp2(1)], 'Location','best')

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);
    FTs1 = FTs(Iv,:);
end

.

[Katie]

Thanks @Star Strider and @Umar

[Star Strider]

My pleasure!

Does it do what you want it to do? I’m still not certain what that is.