Denosing of extremely condensed signal

I am attempting to denoise a high-sampeled spectrum by removing all the sharp peaks and then performing a shape-preserving piecewise cubic spline interpolation to retain as much spectrum information as I can. However, I am not able to clear out all the peaks, including the one peak in the beginning. What am I missing - am I approaching this wisely? Btw., I have several of these spectra that I need to process.

Many thanks for any input in advance!!

Attempt

minval = 0;
maxval = 950;

figure();
plot(lambda, SI, '-')
xlabel('Wavelength')
ylabel('Signal Intensity')
ylim([minval maxval])
title('Original signal')

ind_min = find(SI <= minval);
ind_max = find(SI >= maxval);
lambda([ind_min.', ind_max.']) = [];
SI([ind_min.', ind_max.']) = [];

figure();
plot(lambda, SI, '-')
xlabel('Wavelength')
ylabel('Signal Intensity')
ylim([minval maxval])
title('Signal after removing all values above 950')

% Define threshold representing 25th percentile difference between signal values
diffSI = abs(diff(SI));
thres = quantile(diffSI, 0.25);

%% Remove all sharp peaks
% Remove all signal values differing from neighbouring values above 
% computed threshold and replace their values estimated by 
% shape-preserving piecewise cubic spline interpolation 
for i = 2:(length(SI) - 1)
    if abs(SI(i) - SI(i+1)) > thres
        SI(i) = NaN;
        SI(i+1) = NaN;
    elseif abs(SI(i) - SI(i-1)) > thres
        SI(i) = NaN;
        SI(i-1) = NaN;
    end
end
SI = fillmissing(SI, 'pchip');

figure()
plot(lambda, SI, '-')
xlabel('Wavelength')
ylabel('Signal Intensity')
ylim([minval maxval])
title('Signal after interpolation')

Answer 1

[Image Analyst]

Try movmin() to get envelope. Then try a Saivtky-Golay filter if you want to smooth it out some more.

% Initialization steps:
clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clearvars;
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 15;

s = load('SI_and_lambda.mat')
SI = s.SI;
lambda = s.lambda;

minval = 0;
maxval = 950;
subplot(3, 1, 1);
plot(lambda, SI, '-')
xlabel('Wavelength', 'FontSize', fontSize)
ylabel('Signal Intensity', 'FontSize', fontSize)
ylim([minval maxval])
title('Original signal', 'FontSize', fontSize)
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'Outerposition', [0, 0.05, 1, 0.95]);

% Filter with a moving minimum window.
windowWidth = 27
filteredSignal = movmin(SI, windowWidth);
subplot(3, 1, 2);
plot(lambda, filteredSignal, '-', 'LineWidth', 2)
xlabel('Wavelength', 'FontSize', fontSize)
ylabel('Signal Intensity', 'FontSize', fontSize)
ylim([minval maxval])
title('Filtered signal with moving min', 'FontSize', fontSize)
grid on;

% Filter with a Savitzky-Golay filter.
windowWidthSG = 55
filteredSignal = sgolayfilt(filteredSignal, 2, windowWidthSG);
subplot(3, 1, 3);
plot(lambda, filteredSignal, '-', 'LineWidth', 2)
xlabel('Wavelength', 'FontSize', fontSize)
ylabel('Signal Intensity', 'FontSize', fontSize)
ylim([minval maxval])
title('Filtered signal with moving quadratic Savitzky-Golay filter', 'FontSize', fontSize)
grid on;
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'Outerposition', [0, 0.05, 1, 0.95]);

[De Ar]

Thank you so much for the response! But will the movemin function track the lower envelope as closely as e.g. the imopen filter... ? Because I tried plotting these results on top of each other to see the difference, and the movemin result is slightly downwardly shifted from the spectrum.

[Image Analyst]

movemin() is basically just imerode() - a local min filter.

imclose() is that erosion (local min), followed by a dilation (local max).

I don't know why you'd want a morphological closing since it might introduce more artifacts, but if you feel it gives you a better shape, then go ahead and use it.

The slight downward shift is an artifact that depends on the window size. The bigger you make it, the more the signal is shifted downwards. You might try a modified median filter - like a salt and pepper filter. Basically this filters the signal with a median filter and wherever the signal differs from the median filter by a lot, then it's replaced by the median. See attached demo for a 2-D version. You'd want to adapt it to use medfilt1() instead of medfilt2().