Clear Filters
Clear Filters

Help with torque smoothing and line interpolation

16 views (last 30 days)
Hello, below i have a code snippet fromr a torque smoothing code. Esentially i am looking at torque traces during various conditions. When we apply stimulation, it causes the stim signal to bleed into the torque channel and cause huge spikes. This code snippet graphs the torque area of interest and has the user select the two end points (left and right) from where i want to smooth the data. Esentially you click where the sharp spike occurs, the data between there is deleted and then i have options on how to fill the missing data. The linear intepolation method is ok but obviously leaves a straight line between the data points whereas i was thinking something like a spline interpolation (if that's even the right name for what i want). I have options in the code snippet for higher degree interpolations just in case. But in general i don't really want the inerpolated data to try to follow the original torque signal as it is so far skewed with the bleeding from the stim channel. I hope that makes sense. So in essence i just want the code to look at the two end points and make the transition "more curvy" than just a straight line but at same time i don't want it to follow the original torque trace. You'll see for method 3 in my code that i've commented out several attempts with the help of chat GPT but i'm getting nowhere. Thanks in advance. The specific section i need help with is in Method 3 of the code snippet.
for i = 1:length(VarPulses)
pulse = VarPulses(i);
fig1 = figure;
set(fig1, 'position', get(0, 'ScreenSize'))
% set location of plot titles
ax = axes('position', [0.05, 0.95, 0.92, 0.05]);
xlim([0, time_events(end)])
ylim([0, 5])
text(2, 2.5, figureName, 'FontSize', 16, 'FontWeight', 'bold', 'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle', 'interpreter', 'none')
% Create a string to display the current pulse value
pulseString = ['Current Pulse: ' num2str(pulse)];
text(1, 2.5, pulseString, 'FontSize', 16, 'FontWeight', 'bold', 'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle', 'interpreter', 'none')
set(ax, 'visible', 'off')
% format data plot
ytxt = 'Torque (Nm)';
stim_ytxt = 'Stim (V)';
set(ax, 'XLim', [time_events(1), time_events(end)]);
ybound = get(gca, 'ylim'); % get y-limits of plotted peak
set(zoom, 'Enable', 'on');
% Plot torque events
ax = axes('position', [0.05, 0.1, 0.92, 0.85]);
yyaxis left;
plot(time_events, torq_events(pulse, :))
ylabel(ytxt, 'FontSize', 14)
yyaxis right;
plot(time_events, stim_events(pulse, :))
ylabel(stim_ytxt, 'FontSize', 12)
xlabel('Time(s)', 'FontSize', 14)
set(ax, 'XLim', [time_events(1), time_events(end)]);
% set up variables that will store information about the selected pulses
xtime((i-1)*2 + 1 : i*2, 1) = nan(2, 1);
xidx((i-1)*2 + 1 : i*2, 1) = nan(2, 1);
% user selects two cut points between which to smooth torque
for j = 1:2 % loop twice for two clicks
[x, ~] = ginput(1);
hold on % temporarily holds the current plot and adds a vertical dashed line at the selected x coordinate
line([x x], [ybound(1) ybound(2)], 'linestyle', '--', 'color', [0.9, 0.3, 0.1], 'linewidth', 1.5)
hold off
xtime((i-1)*2 + j, 1) = x;
xidx((i-1)*2 + j, 1) = find(time_events <= x, 1, 'last');
end
%Define start and end indices
start_idx = xidx((i-1)*2 + 1);
end_idx = xidx((i-1)*2 + 2);
% Define the number of points in the cut section
num_points = end_idx - start_idx + 1;
%Cut Torque Interpolations
switch selectedMethod
case {'Method 1: Linear Interpolation'}
% Method 1 = Linear interpolation
interpolated_torque = linspace(torq_events(pulse, start_idx), torq_events(pulse, end_idx), num_points);
case {'Method 2: Higher Degree Polynomial Interpolation'}
% Method 2 = Higher Degree Polynomial Interpolation (2=quadratic, 3=cubic, 4=quartic, 5=quintic...etc.
degree = degreeIdx + 1; % Add 1 to degreeIdx to get the actual degree
x_values = [start_idx, (start_idx + end_idx) / 2, end_idx];
y_values = torq_events(pulse, [start_idx, (start_idx + end_idx) / 2, end_idx]);
x_interp = linspace(start_idx, end_idx, num_points);
coefficients = polyfit(x_values, y_values, degree);
interpolated_torque = polyval(coefficients, x_interp);
case {'Method 3: Spline Interpolation'}
% Method 3 = Spline Interpolation
x_values = start_idx:end_idx;
y_values = torq_events(pulse, start_idx:end_idx);
xi = linspace(start_idx, end_idx, num_points);
% interpolated_torque = interp1(x_values, y_values, xi, 'spline');
% spline_coeffs = spline(x_values, [0, y_values, 0]);
% interpolated_torque = ppval(spline_coeffs,xi);
% interpolated_torque = pchip(x_values, y_values, xi);
alpha = 0.5;
interpolated_torque = (1 - alpha)*y_values(1) + alpha*y_values(end);
end
% Update the interpolated data in the existing variable
torq_events_interp(pulse, start_idx:end_idx) = interpolated_torque;
% Plot torque events with interpolation
ax = axes('position', [0.05, 0.1, 0.92, 0.85]);
yyaxis left;
plot(time_events, torq_events(pulse, :))
hold on
plot(time_events(start_idx:end_idx), interpolated_torque, 'r--', 'LineWidth', 2) % Plots the interpolated data in red dashed line
hold off
ylabel(ytxt, 'FontSize', 14)
yyaxis right;
plot(time_events, stim_events(pulse, :))
ylabel(stim_ytxt, 'FontSize', 12)
xlabel('Time(s)', 'FontSize', 14)
set(ax, 'XLim', [time_events(1), time_events(end)]);
% save current figure (original and smoothed torque)
figurefilename = strcat(PID, exp_ID, 'pulse_', num2str(pulse));
figureFilePath = fullfile(savefilePath, figurefilename);
saveas(fig1, figureFilePath, 'jpg');
close(gcf);
end
  2 Comments
Mathieu NOE
Mathieu NOE on 19 Jan 2024
hello
why not simply replace the noisy signal (around spike location) with a smoothed version of it ?
so simply create the full length smoothed signal and replace the sections you need to
you could use smoothdata , I suggest either with movmedian or gaussian window

Sign in to comment.

Answers (1)

Garmit Pant
Garmit Pant on 30 Jan 2024
Hello Matteo
From what I gather, you are working with torque measurements that are disrupted by sharp spikes from a stimulation signal, and you require a technique to eliminate this noise. The goal is to interpolate the data to maintain the underlying trend and ensure a smooth, natural curve between the user-selected points.
You can refer to the following three methods to achieve the same along with sample code snippets to implement them:
  • Cubic spline interpolation using “fit”: Use the function “fit” and set the “fittype” as “cubicspline”. Use the function “feval” to evaluate the spline over the entire curve.
% Sample data with a spike
x = linspace(0, 10, 100);
y = sin(x); % Adding random noise between x=4 and x=6
y_noise = y + ((x > 4) & (x < 6)) .* (0.5 - rand(1,100));
% Indices of the start and end of the noisy data
start_idx = find(x > 4, 1);
end_idx = find(x < 6, 1, 'last');
% Use data excluding the noisy section to fit the spline
x_fit = [x(1:start_idx-1), x(end_idx+1:end)];
y_fit = [y_noise(1:start_idx-1), y_noise(end_idx+1:end)];
% Fit a cubic spline to the non-noisy data
fit_spline = fit(x_fit(:), y_fit(:), 'cubicspline');
% Evaluate the spline across the entire range, including the noisy section
y_spline = feval(fit_spline, x);
% Plot the results
figure;
plot(x, y_noise, 'b', 'DisplayName', 'Noisy Data');
hold on;
plot(x, y_spline, 'r--', 'DisplayName', 'Cubic Spline Interpolation');
legend;
title('Cubic Spline Interpolation using fit');
xlabel('x');
ylabel('y');
  • PCHIP and Spline Interpolation using “fillmissing”: Replace the data in between the two user selected points with “NaN” values and fill the empty values using “fillmissing”. You can use different fill methods. The code snippet compares “spline” and “phcip” methods.
% Sample data with a spike
x = linspace(0, 10, 100);
y = sin(x); % Adding random noise between x=4 and x=6
y_noise = y + ((x > 4) & (x < 6)) .* (0.5 - rand(1,100));
% Replace the noisy data with NaNs to simulate missing values
y_missing = y_noise;
y_missing((x > 4) & (x < 6)) = NaN;
% Fill missing data using pchip
y_pchip = fillmissing(y_missing, 'pchip', 'SamplePoints', x);
% Fill missing data using spline
y_spline_fill = fillmissing(y_missing, 'spline', 'SamplePoints', x);
% Plot the results
figure;
plot(x, y_noise, 'b', 'DisplayName', 'Noisy Data');
hold on;
plot(x, y_pchip, 'g--', 'DisplayName', 'PCHIP Interpolation');
plot(x, y_spline_fill, 'r--', 'DisplayName', 'Spline Interpolation');
legend;
title('Comparison of PCHIP and Spline Interpolation using fillmissing');
xlabel('x');
ylabel('y');
  • Smoothing the data using “smoothdata”: Use the “smoothdata” function to smooth the spikes between the user selected points. Use “gaussian” smoothing method. The code snippet compares the effects of the “window” input argument. A smaller window size retains details while a larger window size produces smoother data.
% Sample data with a spike
x = linspace(0, 10, 100);
y = sin(x); % Adding random noise between x=4 and x=6
y_noise = y + ((x > 4) & (x < 6)) .* (0.5 - rand(1,100));
smallWindowSize = 5; % Small Gaussian window size
y_smooth_small = smoothdata(y_noise, 'gaussian', smallWindowSize);
% Smooth the data using a large Gaussian window
largeWindowSize = 21; % Large Gaussian window size
y_smooth_large = smoothdata(y_noise, 'gaussian', largeWindowSize);
% Plot the results
figure;
plot(x, y_noise, 'b', 'DisplayName', 'Noisy Data');
hold on;
plot(x, y_smooth_small, 'g--', 'DisplayName', 'Smoothed Data (Small Gaussian)');
plot(x, y_smooth_large, 'm--', 'DisplayName', 'Smoothed Data (Large Gaussian)');
legend;
title('Data Smoothing using smoothdata with Small vs Large Gaussian window');
xlabel('x');
ylabel('y');
For further understanding on the methods mentioned above, you can refer to the following MATLAB Documentation:
  1. “fit” function of Curve Fitting Toolbox: https://www.mathworks.com/help/releases/R2022a/curvefit/fit.html
  2. fillmissing” function: https://www.mathworks.com/help/releases/R2022a/matlab/ref/fillmissing.html
  3. smoothdata” function: https://www.mathworks.com/help/releases/R2022a/matlab/ref/smoothdata.html
I hope you find the above explanation and suggestions useful!

Products


Release

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!