gusucode.com > 信号处理工具箱 - signal源码程序 > signal\signal\signal\sgolayfilt.m
function y=sgolayfilt(x,k,F,varargin) %SGOLAYFILT Savitzky-Golay Filtering. % SGOLAYFILT(X,K,F) smoothes the signal X using a Savitzky-Golay % (polynomial) smoothing filter. The polynomial order, K, must % be less than the frame size, F, and F must be odd. 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 K equals F-1, no smoothing % will occur. % % SGOLAYFILT(X,K,F,W) specifies a weighting vector W with length F % containing real, positive valued weights employed during the % least-squares minimization. % % See also SGOLAY, MEDFILT1, FILTER % References: % [1] Sophocles J. Orfanidis, INTRODUCTION TO SIGNAL PROCESSING, % Prentice-Hall, 1995, Chapter 8. % Author(s): R. Losada % Copyright (c) 1988-98 by The MathWorks, Inc. % $Revision: 1.3 $ $Date: 1998/07/21 11:53:16 $ error(nargchk(3,4,nargin)); % Check if the input arguments are valid if round(F) ~= F, error('Frame length must be an integer.'), end if rem(F,2) ~= 1, error('Frame length must be odd.'), end if round(k) ~= k, error('Polynomial degree must be an integer.'), end if k > F-1, error('The degree must be less than the frame length.'), end [mx,nx]=size(x); if mx == 1, x = x(:); nrows = nx; ncols = 1; else nrows = mx; ncols = nx; end if nrows < F, error('The length of the input must be >= frame length.'), end if nargin < 4, % No weighting matrix, make W an identity W = ones(F,1); else W = varargin{1}; % Check for right length of W if length(W) ~= F, error('The weight vector must be have the same length as the frame length.'),end % Check to see if all elements are positive if min(W) <= 0, error('All the elements of the weight vector must be greater than zero.'), end end % Compute the projection matrix B B = sgolay(k,F,W); % Compute smoothing result for each signal column: for col=1:ncols, % Compute the transient on wit = flipud(x(1:F,col)); yit = flipud(B((F-1)./2+1:end,:))*wit; % Compute the steady state output [xb,z] = buffer(x(:,col), F, F-1, 'nodelay'); % Buffer the input xb = xb(:,2:end-1); % Take the part used for steady state yss = (B((F-1)./2+1,:)*xb).'; % Compute the steady state output % Compute the transient off wot = flipud(x(end-(F-1):end,col)); yot = flipud(B(1:(F-1)./2+1,:))*wot; % Form the total output y(:,col) = [yit;yss;yot]; end