Denosing of extremely condensed signal

2 次查看（过去 30 天）
De Ar
De Ar 2018-12-7
评论： Image Analyst 2018-12-18
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!!
SI_original.png
SI_after_950.png
SI_after_interp.png
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')

请先登录，再回答此问题。

采纳的回答

Image Analyst
Image Analyst 2018-12-7
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]);
0000 Screenshot.png
  2 个评论
De Ar
De Ar 2018-12-18
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
Image Analyst 2018-12-18
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().

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Signal Generation and Preprocessing 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by