Savitzky-Golay smooth/differentiation filters and filter application

Routine to generate Savitzky-Golay smoothing and differentiation filters and routine to apply these

Description

Contains savitzkyGolay.m and savitzkyGolayFilt.m:
Function:
        Savitzky-Golay Smoothing and Differentiation Filter
        The Savitzky-Golay smoothing/differentiation filter (i.e., the
        polynomial smoothing/differentiation filter, or the least-squares
        smoothing/differentiation filters) optimally fit a set of data
        points to polynomials of different degrees.
        See for details in Matlab Documents (help sgolay). The sgolay
        function in Matlab can deal with only symmetrical and uniformly
        spaced data of even number.
        This function presented here is a general implement of the sgolay
        function in Matlab. The Savitzky-Golay filter coefficients for even
        number, nonsymmetrical and nonuniformly spaced data can be
        obtained. And the filter coefficients for the initial point or the
        end point can be obtained too. In addition, either numerical
        results or symbolical results can be obtained. Lastly, this
        function is faster than MATLAB's sgolay.
 
  Usage:
        [fc,df] = savitzkyGolay(x,n,dn,x0,flag)
    input:
        x = the original data point, e.g., -5:5
        n = polynomial order
        dn = differentation order (0=smoothing), default=0
        x0 = estimation point, can be a vector default=0
        W = weight vector, can be empty
               must have same length as x0 default=identity
        flag = numerical(0) or symbolical(1), default=0
 
    output:
        fc = filter coefficients obtained (B output of sgolay).
        df = differentiation filters (G output of sgolay).

.

.

savitzkyGolayFilt(X,N,DN,F) filters the signal X using a Savitzky-Golay
    (polynomial) filter. The polynomial order, N, must be less than the
    frame size, F, and F must be odd. DN specifies the differentiation
    order (DN=0 is smoothing). For a DN higher than zero, you'll have to
    scale the output by 1/T^DN to acquire the DNth smoothed derivative of
    input X, where T is the sampling interval. The length of the input X
    must be >= F. If X is a matrix, the filtering is done on the columns
    of X.
 
    Note that if the polynomial order N equals F-1, no smoothing
    will occur.
 
    savitzkyGolayFilt(X,N,DN,F,W) specifies a weighting vector W with
    length F containing real, positive valued weights employed during the
    least-squares minimization. If not specified, or if specified as
    empty, W defaults to an identity matrix.
 
    savitzkyGolayFilt(X,N,DN,F,[],DIM) or savitzkyGolayFilt(X,N,DN,F,W,DIM)
    operates along the dimension DIM.

Acknowledgements

The author wishes to acknowledge the following in the creation of this submission:
Savitzky-Golay Smoothing and Differentiation Filter
This submission has inspired the following:
iSignal 2.0

• MATLAB release: MATLAB 7.11 (2010b)

[Omar]

Just a question, but isn't the output (y) "mirrorred" when you take the first derivative of the signal, i.e y=savitzkyGolayFilt(X,2,1,11); X being the unfiltered signal?

[Diederick]

It is, its sign is opposite. This is consistent with Matlab's sgolay function

[Brandon Teska]

It's an easy fix, I just add the following to the bottom of this module:

% Flip data for odd derivatives
if mod(DN,2)
    y = -y;
end

Thanks for writing this, it's a great m-file.

[Diederick]

Thank you Brandon :)
That fix would do the trick. But I prefer staying consistent with matlab's sgolay so this can be used as a more advanced drop-in replacement. Thanks for documenting it here though, as others might run into the same issues!

[rjones]

 

[David McElroy]

G'Day,

I'm using a 4th derivative with a sampling interval of 26. Could we see an example of how to " scale the output by 1/T^DN"?

Thanks!

[Diederick]

Hi David,

I assume you mean that your data is sampled at 26 Hz?

If you get your filter output as follows
output = savitzkyGolayFilt(input,polynomial_order,4,filter_length);

then you would have to scale as follows:
output_correct_scaling = output*26^4;
which is equivalent to:
T = 1/26; % sampling interval: 1/26th second between each sample
output_correct_scaling = output/T^4;

Hope that helps,
Best,
Dee

[David McElroy]

Thanks for the quick reply, Dee. My data is wavelength reflectance data which is equally spaced from 350-2500nm, so my sampling interval is in nm. On running your suggestion, it appears the scaling itself was apparently not the issue - the data points are sampled 1nm apart, but I want to smooth at the level of 26nm. Is that possible with some call to this function, or what would I need to use something else?

[Diederick]

Hi David,

I'm actually not sure what you mean by "smooth at the level of 26nm". Maybe its better to ask on the matlab mailing list (see newsgroup in the navigation menu).

Best,
Dee