eigen_filtering.m

著名的seismiclab的代码 是地震学研究人员必备的工具

 function [s,w,g] = eigen_filtering(y,p,mu);

% Given a 1D noisy sequence y, the order p of
% the ARMA(p,p) model and the regularization parameter mu
% this function computes the clean signal s, an estimate of 
% the noise sequence, and the predicor error fiter g.

 
 N = length(y);                     
 Y = convmtx(y,p+1);                 % data convolution matrix
 R = Y'*Y/N;                         % data correlation matrix
 [g,Pw] = eigs(R,1,'SM');            % compute 1 eigenvalue (the SMallest)
 g0 = g(1)
 g  = g/g0;
 e = Y*g;                            % non-white sequence e
 G = convmtx(g,N);                   % correlation matrix of the filter
 D = eye(N)*mu;                      % regularization term
 w = inv(G'*G+D)*G'*e;               % estimate of the noise
 s = y-w;                            % estimate of the clean signal

 return