How can I apply sgolay derivatives to a matrix of signals

3 次查看(过去 30 天)
Dear all,
My first post here. Thank you for helping me out!
I am trying to calculate derivatives of a signal, using sgolay filtering, as the signals are slightly noisy. I have a working function, see below.
However, this uses loops to loop through the individual signals. I would like to create it as a matrix operation, if this makes sense from a performance point of view. The signals are typically 2000 elements long. A dataset may contain a few hundred of these signals. I will use the script on different places so it may have 20 implementations in the final workflow.
So, the questions are:
  • Does it make sense to rework this to a matrix operation?
  • Can someone help me with the re-work of the code? I am absolutely lost..
Current script:
function [derivatives] = JF_savitzky_derivative(spectra, order, framelen, returnorder)
% create return data array
derivatives = spectra;
derivatives(:,:) = 0;
% define input
sizes = size(spectra);
numspectra = sizes(1);
numbands = sizes(2);
% Create the filter
[b,g] = sgolay(order,framelen);
% create the temporary valuestore
dx = zeros(numbands,returnorder+1);
% Set the spacing between measurements
dt = 1;
% for each spectrum run the filter
for spectrum = 1:length(spectra(:,1))
x=spectra(spectrum,:);
for p = 0:returnorder
dx(:,p+1) = conv(x, factorial(p)/(-dt)^p * g(:,p+1), 'same');
end
derivatives(spectrum,:) = dx(:,returnorder+1);
end
end

采纳的回答

Grant Sellers
Grant Sellers 2018-1-25
Hello Jelle,
I think it does make sense to rework your script as a matrix operation. Typically, vectorization does have good performance improvements when working with larger arrays like this. Even if performance is not a desire, it often makes the code shorter and easier to read.
As far as how to do that, it seems like you will want to use the "conv2" function to convolve your matrix of spectra with the filter vector. This will get rid of the for loop that iterates over spectra. I think what you will want to do next is to replace the second for loop that iterates over the order of the filter response with something simpler like p=return order. You are not currently using any of the dx terms except the last, and since they are not used in computation, they will not need to be calculated.
If you do need each order of the derivative for each signal, then I think you will need to keep one for loop that iterates over p that calls the "conv2" function on the spectra matrix and the filter vector.
One more note, you mentioned wanting performance improvements with this code. Consider replacing
% create return data array
derivatives = spectra;
derivatives(:,:) = 0;
with
derivatives=zeros(size(spectra));
This implementation will save time as spectra grows larger.
  2 个评论
Jelle
Jelle 2018-1-29
Hi Grant,
This seems like a way to go. As you may have noticed, I am very new at matlab so trying to get my head around all this stuff is really time consuming. Very glad that there are persons out there like yourself, willing to help newbies on there way! Thank you!
I will try to get this implemented maybe today else tomorrow and let you know how I go.
Jelle
Jelle 2018-1-29
Brilliant!
% Transpose the data array
x=spectra';
% Create the filter
[b,g] = sgolay(order,framelen);
% Set the spacing between measurements
dt = 1;
% 2d convolution filter to get the derivative at the desired order
derivatives = conv2(x, factorial(returnorder)/(-dt)^returnorder * g(:,returnorder+1), 'same');
all works & gives me the output I expect. :)
Thank you!

请先登录,再进行评论。

更多回答(0 个)

Community Treasure Hunt

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

Start Hunting!

Translated by