Help with torque smoothing and line interpolation

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

[Mathieu NOE]

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

[Mathieu NOE]

hello again

Problem solved ?

Answer 1

[Garmit Pant]

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!