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.

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.

.