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