How can I apply sgolay derivatives to a matrix of signals

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

Answer 1

[Grant Sellers]

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.

[Jelle]

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]

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!