Determining optimal cut off frequency

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.

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.

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.

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.