golub1.m

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

function [theta, X1] = golub1(X1, M1, P, s)
 
%%MATLAB code to estimate the p eigenvectors associated with the p
%%largest eigenvalues of the n-by-n matrix A = P*M*P, where:
%%input:   X1 =  orthogonal [n-by-p] matrix; original estimate of
%%               sought eigenvectors
%%         M1 = first row of the Toeplitz matrix M
%%         P  = [n-by-n] projection matrix: A=P*M*P
%%         s  = block dimension of the Lanczos procedure
%%output:  theta = estimated p largest eigenvalues of A
%%         X1    = associated p eigenvectors of A
%%example1: prolate tapers with bandwidth w=0.0625
%%          N=64; p=5; I=ones(N, 1); w=0.0625;
%%          x=gst(N, p, I); 
%%          temp=(0:(N-1)); M1=sin(2*pi*w*temp)./(pi*temp); M1(1)=w;
%%          P=diag(I); s=3;
%%          [theta, prolx]=golub1(x, M1, P, s);
%%example2: minimum bias tapers
%%          N=64; p=5; I=ones(N, 1); 
%%          x=gst(N, p, I); P=diag(I); s=3;
%%          M1=((-1).^(1:N))./(2*pi*pi*(0:(N-1)).^2); M1(1)=11/12;
%%          [theta, mbx]=golub1(x, M1, P, s);
 
p=length(X1(1,:));
count=1;
while (count > 0)
        X1=lanczos1(X1, M1, P, s);
        Tbar=P*X1;
        Tbar=cycmult(M1, Tbar);
        Tbar=P*Tbar;
        Tbar = X1'*Tbar;
        [S,T]=schur(Tbar);
        [theta, in] = mysort(diag(T));
        S = S(:, in);
        T = diag(theta);
        theta=theta(1:p); 
 
        Res=P*X1;
        Res=cycmult(M1, Res);
        Res=P*Res-X1*T;
        Res=Res-X1*(X1'*Res);
        Res=Res(:, (1:p));
        X1 = X1*S(:, (1:p));
 
        count=0;
        for i=1:p
                if (norm(Res(:,i)) > eps)
                        count=count+1;
                end
        end
 
end