To obtain the cumulative signal power spectrum from the single-sided spectrum of a biomechanical signal to identify the frequency under which 99% of the signal power is

I am looking to select an appropriate cut-off frequency for a Butterworth filter applied to biomechanical data. I have read the process of performing a fast fourier transformation to take the signal from a time to frequency domain before then identifying the frequency under which 99% of the signal power is contained to quantify the cutoff frequency. Below is an image I have frequently seen to outline this process.

I have come across some resources to help with this process but it is a bit of a minefield when you are attempting to learn MATLAB at the same time as FFT. Has anyone got any advice or tips towards completing this in MATLAB? If I can plot the cumulative signal power % against the frequency then I can identify the appropriate cutoff frequency for a Butterworth filter as the frequency for which 99% of the signal power resides under.

Many thanks,

T

[Alexander]

Honestly, I don't understand your problem. What are the time data and what is the samplerate?

[Thomas Eli]

I have ground reaction force data collected from force platforms sampling at 1000Hz over a 10 second period.

[Star Strider]

That Fourier transform doesn’t look correct to me. It must have been processed in some (currently opaque) way.

I use the envelope function to ‘outline’ the fft results, and then use the upper envelope result to design filters.

It would definitely help to have more information, especially the signal, a description of the analysis approach (that is not obvious to me), and the code currently used to analyse the data.

[Alexander]

If you are sampling at 1 kHz the input filter must be lower than 500 Hz cutoff frequency (theoretical) and 100-200 (engineering approach). Why don't you make a simple fft of your data? Mirrored frequencies should be clearly visible, especially if do an oversampling with (e.g.) 10 kHz. In the most cases (my experience, automotive) there were no filters necessary. What kind of excitation do you have?

[Thomas Eli]

Below is the code currently used and signal attached - all I have it doing thus far is essentially importing the data of interest which is from two force platforms that are summed to create 'Fz' (the signal of interest).

The following generally tends to be the information provided as to how/what is done - "The first option is to perform at spectral analysis of the marker trajectories using a fast fourier transform (FFT) to examine the cumulative content of the signal in the frequency domain (Figure 1 - this is the Figure I attached in my original post). Typically the choice of cut-off is taken as the frequency at which either 95 or 99% of the signal power is contained below.

Below is the paper the above is from...

Sinclair, J., Taylor, P. J., & Hobbs, S. J. (2013). Digital filtering of three-dimensional lower extremity kinematics: An assessment. Journal of human kinetics, 39, 25.

%% Clear Workspace & Set Directory
clear; % clear 'leftover' variables from the workspace
clc; % clear previous text or output from the command window
close all; % close all windows or figures that were open

% Set working folder & file
myDir = uigetdir; % open folder to select working directory
myFiles = dir(fullfile(myDir,'AEL_CMJ_1123.txt')); % locating CMJ text file - change every analyses

% Import CMJ trial
BaseFileName = myFiles(1).name; 
FullFileName = fullfile(myDir, BaseFileName); % builds a full filename
fprintf(1,'Now reading %s\n',FullFileName); % display that the file is being read by the programme
num = readmatrix(FullFileName); % import the file

% Extract the data
FzR = num(:,1); % Fz data from force plate 1 (right leg)
FzL = num(:,2); % Fz data from force plate 2 (left leg)
Fz = FzR+FzL; % Fz data from both force plates summed

Answer 1

[Star Strider]

I cannot find a time vector or sampling frequency. I need at least one of those.

The article is open-source, however not illuminating with respect to the signal processing.

This should get you started —

M1 = readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1339494/AEL_CMJ_1123.txt')

M1 = 10000×2

387.3942 387.6488 388.0743 386.6169 386.5396 388.6821 388.0718 386.7832 387.3948 387.3062 387.3934 386.6183 387.7338 386.6205 388.4194 386.6141 387.3925 386.2696 388.2491 386.6045

Fs = 1000;                                      % Sampling Frequency (Default Value)
L = size(M1,1);
t = linspace(0, L-1, L).'/Fs;                   % Time Vector

figure
plot(t, M1)
grid
xlabel('Time (s)')
ylabel('$Amplitude (\frac{g\cdot cm}{s^2})$', 'Interpreter','latex')
legend('R','L', 'Location','best')

Fn = Fs/2;
NFFT = 2^nextpow2(L)

NFFT = 16384

FTM1 = fft((M1-mean(M1)).*hann(L), NFFT)/L;
Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);

figure
plot(Fv, abs(FTM1(Iv,:))*2)
grid
xlabel('Frequuency')
ylabel('Magnitude')
legend('R','L', 'Location','best')
xlim([0 0.2*Fs])

ctz = cumtrapz(Fv,(abs(FTM1(Iv,:))*2).^2);
figure
plot(Fv, ctz)
grid
xlabel('Frequuency')
ylabel('Cumulative Power')
legend('R','L', 'Location','best')
xlim([0 0.2*Fs])

ctzn = ctz./ctz(end,:);

cpnv = [0.9; 0.95; 0.99];
for k = 1:size(ctzn,2)
    Fv_vals(:,k) = interp1(ctzn(:,k), Fv, cpnv);
end
% Fv_vals

figure
plot(Fv, ctzn)
grid
xlabel('Frequuency')
ylabel('Cumulative Power (Normalised)')
legend('R','L', 'Location','best')
xlim([0 0.2*Fs])

FrequencyLimits = table(cpnv*100,Fv_vals(:,1),Fv_vals(:,2), 'VariableNames',{'PowerVal (%)','Right (Hz)','Left (Hz)'})

FrequencyLimits = 3×3 table
    PowerVal (%)    Right (Hz)    Left (Hz)
    ____________    __________    _________

         90           6.3456       6.2923  
         95           9.8479        9.281  
         99           36.201       41.388  

EDIT — (29 Mar 2023 at 15:12)

Updated sampling frequency, re-ran code

.

[Thomas Eli]

Thank you for sending this across. The sampling frequency of the data in this file/analysis is 1000Hz - does this help?

[Star Strider]

My pleasure!

It does. I updated the sampling frequency and then re-ran my code.

I am still not exactly certain what you want to do. An acceptable baandpass frequency design would be:

Fz_Filt = lowpass(M1, mean(FrequencyLimits{end,[2 3]}), Fs, 'ImpulseResponse','iir');

And here that produces:

M1 = readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1339494/AEL_CMJ_1123.txt')

M1 = 10000×2

387.3942 387.6488 388.0743 386.6169 386.5396 388.6821 388.0718 386.7832 387.3948 387.3062 387.3934 386.6183 387.7338 386.6205 388.4194 386.6141 387.3925 386.2696 388.2491 386.6045

Fs = 1000;                                      % Sampling Frequency (Default Value)
L = size(M1,1);
t = linspace(0, L-1, L).'/Fs;                   % Time Vector

figure
plot(t, M1)
grid
xlabel('Time (s)')
ylabel('Amplitude $(\frac{g\cdot cm}{s^2})$', 'Interpreter','latex')
legend('R','L', 'Location','best')

Fn = Fs/2;
NFFT = 2^nextpow2(L)

NFFT = 16384

FTM1 = fft((M1-mean(M1)).*hann(L), NFFT)/L;
Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);

figure
plot(Fv, abs(FTM1(Iv,:))*2)
grid
xlabel('Frequency')
ylabel('Magnitude')
legend('R','L', 'Location','best')
xlim([0 0.2*Fs])

ctz = cumtrapz(Fv,(abs(FTM1(Iv,:))*2).^2);
figure
plot(Fv, ctz)
grid
xlabel('Frequency')
ylabel('Cumulative Power')
legend('R','L', 'Location','best')
xlim([0 0.2*Fs])

ctzn = ctz./ctz(end,:);

cpnv = [0.9; 0.95; 0.99];
for k = 1:size(ctzn,2)
    Fv_vals(:,k) = interp1(ctzn(:,k), Fv, cpnv);
end
% Fv_vals

figure
plot(Fv, ctzn)
grid
xlabel('Frequency')
ylabel('Cumulative Power (Normalised)')
legend('R','L', 'Location','best')
xlim([0 0.2*Fs])

FrequencyLimits = table(cpnv*100,Fv_vals(:,1),Fv_vals(:,2), 'VariableNames',{'PowerVal (%)','Right (Hz)','Left (Hz)'})

FrequencyLimits = 3×3 table
    PowerVal (%)    Right (Hz)    Left (Hz)
    ____________    __________    _________

         90           6.3456       6.2923  
         95           9.8479        9.281  
         99           36.201       41.388  

Fz_Filt = lowpass(M1, mean(FrequencyLimits{end,[2 3]}), Fs, 'ImpulseResponse','iir');

figure
plot(t, Fz_Filt)
grid
xlabel('Time (s)')
ylabel('Amplitude $(\frac{g\cdot cm}{s^2})$', 'Interpreter','latex')
legend('R','L', 'Location','best')

The filter doesn’t do much filtering because of the relatively high cutoff frequency. By design, it lets 99% of the signal power through.

.

[Alexander]

I'm a little bit late(it took some time to analyse), and I think that the most things are mentioned in the answer above, I'll post it anyway so the work was not for nothing.

My short analyse of your data, w/o nowing your test setup, test and equipment. I assume the dimension of your data is N. It seems there is a base load of about 400 N. After 3 s the load decreaes to 800N and then something brakes(?). After that a very fast drop to (low) negative forces, as something is pulling at the force transducer. This remains for half a second and then raising again very quickly. Sometimes FzL and FzR are showing a diametral behavior (good seen after 5s). Wobbling platform?

My steps:

1. Fz I set to (FzR+FzL)/2. I think a simple sum makes no sense, but that is your decision. I'm now focused on FZ.

2. Looking at the time data in figure(1) what I have assumed above.

3. In figure(2) (below) you can see the pure fft of Fz (not much).

4. In figure(3) I removed the mean of Fz to make the interesting part more visible. The mean is allways around the frequency of 0Hz and scales the y-axis.

Some remarks to "my" fft I used. I don't like normalized data and also a amplitude density in 1/Hz is not very interpretable on the first view. In the plots attached you see frequency in Hz and the amplitudes not amplitudes/Hz. This helps to get a quick first view, but you can't use this for a ifft. I will attach the data so you can make your plots yourselfe as you like.

IMHO, you don't need a filter. Of cause this depends on the noise and frequencies you want to accept. You will find in the plots and data, above 250Hz that more or less nothing happens.

If you want to have a filter anyway, you can use a butterworth whith 100Hz cutoff frequency and order of 2 (sample frequency 1 kHz) with the coefficients:

B =

0.0675 0.1349 0.0675

A =

1.0000 -1.1430 0.4128

Frequency response and phase in figure(4).

If you have any questions or you find a mistake: please report.