sgsmooth.m

这个帖子中我想讨论的是移动窗口多项式最小二乘拟和平滑方法

function y=sgsmooth(x,width,order)

%   A smooth method first mentioned by Savitzky and Golay in the paper
%   Abraham Savitzky & Marcel J. E. Golay. Smoothing and Differentiation of Data by 
%   Simplified Least Squares Procedures. Analytical Chemistry.1964,36(8):1627-1639.
%   
%   Then corrected by Jeans Teinier, Yves Termonia & Jules Deltour in 1972,
%   published by Analytical Chemistry as well. And Peter A. Gorry in his
%   article: General Least-Squares Smoothing and Differentiation by the Convolution 
%            (Savitzky-Golay) Method (Anal.Chem.1990,62(6): 570-573) showed us the
%   generalized implementation of SG least square smoothing method.
%
%   The algorithm showed in this program mainly from Jeans Ternier etc.for calculating 
%   the weighting parameters ,also from Peter A. Gorry for calculating the leading and 
%   trailing value.Here we set the derivative order s=0. 
%   
%   y=SavGo(x,width,order):
%   x ---------- Data you must input and it must be a vector.
%
%   width ------ The window width for smoothing.
%
%   order ------ The order of polynomials. It must be smaller than the window width you 
%                choose for smooothing.
%
%   DNP. 2007.12.16

y=[];
[m,n]=size(x);

if ~any([m,n]==1)
    error('The data you want to be smoothed should be a vector!');
    return
elseif m==1
    x=x';              % Make x be a column vector.
end

if width<=order
    error('The order of polynomials must be smaller then the window width you choose!');
    return
elseif all([width order]<=0)
    error('The window width and order must be positive!');
    return
end

if nargin < 3
    error('You must input data for smoothing, and choose the window width and polynomial order.');
    return
elseif nargin > 3
    error('Too many inputs!');
    return
end

if mod(width,2)==0
    error('The window width should be odd.');
    return
end

m=(width-1)/2;
B=[];
h=zeros(width);

% Coefficient matrix ................
for i=-m:m
%     for t=-m:m
        for o=0:order
%             coef=((2*o+1)*factorial(2*m)/factorial(2*m-o))/(factorial(2*m+o+1)/factorial(2*m));
%             p_i=grampoly(i,o,m);
%             p_t=grampoly(t,o,m);
%             h(i+m+1,t+m+1)=h(i+m+1,t+m+1)+coef*p_i*p_t;
            B(i+m+1,o+1)=i^o;
%         end
    end
end

% ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
% ......... Smoothing data ..............

% y=zeros(length(x),1);
for i=1:length(x)
    if i-1-m < 0
        y(i)=h(:,i)'*x(1:2*m+1);
    elseif i > length(x)-m
        y(i)=B(2*m+1-(length(x)-i),:)*(br\bq')*x(length(x)-2*m:end);
    else
        [bq br]=qr(B,0);         % QR decomposition of B for Least Square Fitting.
        y(i)=B(m+1,:)*(br\bq')*x(i-m:i+m);
    end
end
        
        
                                                                                                                                                
                
        
    
% % **************************************************************************************
% % ............ Gram polynomials. .............
% 
% function p=grampoly(t,k,m)
% 
% %   t ------- The position for calculation.
% %   k ------- The polynomial order.
% %   m ------- The half of window width (2m+1).
% 
% p=[];
% 
% if k==0
%     p=1;
% elseif k>0
%     p(1)=0;   % p(0)=1, p(-1)=0
%     p(2)=1;
%     for j=1:k
%         p(j+2)=(2*(2*j-1)/(j*(2*m-j+1)))*t*p(j)-((j-1)*(2*m+j)/(j*(2*m-j+1)))*p(j+1);
%     end
% end
% 
% p=p(end);