均匀采样信号的鲁棒Savistky-Golay滤波(MATLAB）

S-G滤波器又称S-G卷积平滑器，它是一种特殊的低通滤波器，用来平滑噪声数据。该滤波器被广泛地运用于信号去噪，采用在时域内基于多项式最小二乘法及窗口移动实现最佳拟合的方法。与通常的滤波器要经过时域－频域－时域变换不同，S-G滤波直接处理时域数据进行平滑，其平滑效果随窗口宽度不同而不同。相对于均值平滑滤波，S-G滤波更能保留相对极大值、极小值和宽度等分布特征。该滤波算法的另一优点是其运算量相对较小，对计算机的内存及数据处理能力要求较低。鉴于此，采用一种鲁棒Savistky-Golay滤波器对均匀采样信号进行降噪平滑，运行环境MATLAB 2021。

for i = 1:step:Nind1 = max(1,i-half_window) ;ind2 = min(N,i+half_window) ;y_spl = y_input(ind1:ind2) ;x_spl = [ind1:ind2]' ;weights = ones(size(y_spl)) ;p = inf(order+1,1) ; s = inf ; ds = 1 ; dp = 1 ;X = repmat(x_spl,1,order+1).^[order:-1:0] ;c = 4.685 ; nb_iter = 0 ;while (abs(ds) >= opts.TolFun) || (abs(dp) >= opts.TolX)if nb_iter > opts.MaxIterbreakendnb_iter = nb_iter+1 ;old_p = p ;p = (sqrt(weights).*X)\(sqrt(weights).*y_spl) ;y_spl_calc = X*p ;r = y_spl_calc-y_spl ;tau = median(abs(r-median(r)))/0.6745 ;if tau == 0s = std(r)*sqrt(1+1/numel(y_spl)+(i-mean(x_spl)).^2/sum((x_spl-mean(x_spl)).^2)) ;breakendz = r/tau ;weights = (abs(z)<c).*(1-(z/c).^2).^2 ;weights = weights.*(weights>=0) ; old_s = s ;s = std(r,weights)*sqrt(1+1/numel(y_spl)+(i-sum(weights.*x_spl)./sum(weights)).^2/sum((x_spl-sum(weights.*x_spl)./sum(weights)).^2)) ;dp = max((old_p-p)./p) ;ds = (old_s-s)/s ;endy3(i,:) = p' ;y2(i,1) = s ;y1(i,1) = polyval(p,i) ;endif y_input_type == "row"y1 = y1' ;y2 = y2' ;x3 = y3' ;y_input = y_input' ;x_interp = 1:step:N ;x_query = 1:N ;
elsex_interp = [1:step:N]' ;x_query = [1:N]' ;
endif step == 1% No need for interpolation.return
endy1 = interp1(x_interp,y1(1:step:N),x_query,opts.Interpolation,"extrap") ;
y2 = interp1(x_interp,y2(1:step:N),x_query,opts.Interpolation,"extrap") ;
y3 = interp1(x_interp,y3(1:step:N,:),x_query,opts.Interpolation,"extrap") ;

• 工学博士，担任《Mechanical System and Signal Processing》《中国电机工程学报》《控制与决策》等期刊审稿专家，擅长领域：现代信号处理，机器学习，深度学习，数字孪生，时间序列分析，设备缺陷检测、设备异常检测、设备智能故障诊断与健康管理PHM等。代码通过知乎学术咨询获得：https://www.zhihu.com/consult/people/792359672131756032?isMe=1