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.
