gst.m

工程计算MATLAB code to calculate the reorthogonalized sine tapers input: N = the length of the time se

function [X] = gst(N, p, I)

%%MATLAB code to calculate the reorthogonalized sine tapers
%%input:   N = the length of the time series data to be tapered
%%         p = the number of tapers requested
%%         I = the gap structure; a vector of length N 
%%            I(t) = 1 if there is data at time t, t=1, ..., N
%%            I(t) = 0 if there is a gap at time t
%%output:  X = N-by-p vector of the reorthogonalized sine tapers      
%%example: N=64; p=5; I=ones(N, 1); 
%%         x=gst(N, p, I); 

B=(1:N)';
V=spdiags(B,0,N,N)*(1+sparse(N,N));
B=(1:p)';
Q=V*[spdiags(B,0,p,p);sparse(N-p, p)];
Q=Q*pi/(N+1);
Q=sqrt(2/(N+1))*sin(Q);
 
P=spdiags(I,0,N,N); 
X=P*Q;
X=X*(X'*X)^(-.5);