savgol.m

偏最小二乘算法在MATLAB中的实现

function y_hat = savgol(y,width,order,deriv)
%SAVGOL Savitsky-Golay smoothing and differentiation
% Inputs are the matrix of row vectors to be smoothed (y),
% and the optional variables specifying the number of points in
% filter (width), the order of the polynomial (order), and the
% derivative (deriv). The output is the matrix of smoothed
% and differentiated row vectors (y_hat). If number of points,
% polynomial order and derivative are not specified,
% they are set to 15, 2 and 0, respectively. I/O format is:
% y_hat = savgol(y,width,order,deriv);
%
% Example: if y is a 5 by 100 matrix
%          then savgol(y,11,3,1) gives the 5 by 100 matrix of
%          first-derivative row vectors resulting from a 11-
%          point cubic Savitzky-Golay smooth of each row of y

% Sijmen de Jong
% Unilever Research Laboratorium Vlaardingen
% Feb 1993
% Modified by Barry M. Wise
% May 1994
[m,n] = size(y);
y_hat = y;
% set default values: 15-point quadratic smooth
if nargin<4
  deriv= 0;
  disp('  '), disp('Derivative set to zero') 
end
if nargin<3
  order= 2; 
  disp('  '), disp('Polynomial order set to 2')
end
if nargin<2
  width=min(15,floor(n/2)); 
  s = sprintf('Width set to %g',width);
  disp('  '), disp(s)  
end
% In case of input error(s) set to reasonable values
w = max( 3, 1+2*round((width-1)/2) );
if w ~= width
  s = sprintf('Width changed to %g',w);
  disp('  '), disp('Width musth be >= 3 and odd'), disp(s)
end
o = min([max(0,round(order)),5,w-1]);
if o ~= order
  s = sprintf('Order changed to %g',o); disp('  ')
  disp('Order must be <= width -1 and <= 5'), disp(s)
end
d = min(max(0,round(deriv)),o);
if d ~= deriv
  s = sprintf('Derivative changed to %g',d); disp('  ')
  disp('Deriviative must be <= order'), disp(s)
end
p = (w-1)/2;
% Calculate design matrix and pseudo inverse
x = ((-p:p)'*ones(1,1+o)).^(ones(size(1:w))'*(0:o));
for k = 1:m
  weights = (x'*x)\x';                          
% Smoothing and derivative for bulk of the data
  for i=p+2:n-p-1                              
    y_hat(k,i) = weights(d+1,:)*y(k,i-p:i+p)';      
  end
  % Smoothing and derivative for tails 
  weights = weights*[y(k,1:w)', y(k,n-w+(1:w))']; % full polynomial model
  for i=1:d
    weights = diag(1:o+1-i)*weights(2:o+2-i,:); % or its d'th derivative
  end
  y_hat(k,1:p+1)     = (x(1:p+1,1:1+o-d)*weights(:,1))'; % fitting the tails
  y_hat(k,n-p+(0:p)) = (x(p+1:w,1:1+o-d)*weights(:,2))';
end