How can I apply sgolay derivatives to a matrix of signals

Illustration
technicalsouce9 - 2021-01-14T15:26:25+00:00
Question: How can I apply sgolay derivatives to a matrix of signals

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

Expert Answer

Profile picture of Prashant Kumar Prashant Kumar answered . 2025-11-20

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.


Not satisfied with the answer ?? ASK NOW

Get a Free Consultation or a Sample Assignment Review!