readme

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

*********************************************************	
*							*
*	MATLAB FUNCTIONS				*
*							*
*********************************************************



gst.m
-----
	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); 


golub1.m
--------
	
	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);
 

subroutines needed in the golub1.m function
-------------------------------------------

lanczos1.m
----------
	
	MATLAB code to implement the block Lanczos procedure
	input:  X1 = [n-by-p] initial guess for p eigenvectors of A=P*M*P
	        M1 = [1-by-n] the first row of the Toeplitz matrix M
	        P  = [n-by-n] projection matrix P
		s  = block dimension of the Lanczos procedure
	output: X1 = [n-by-(s+1)*p] matrix obtained by the block Lanczos
	             procedure 
 

cycmult.m
--------- 
	
	MATLAB function  
	input:   a = vector of length n
	         X = [n-by-p] matrix
	output:  R = [n-by-p] = M*X
	             where M is [n-by-n] symmetric Toeplitz matrix whose
	             first row is the input vector a
 

mysort.m
---------

	MATLAB function 
	input:  x  = vector 
	output: y  = sorted in decreasing order x
	        in = the index of the original x values; x(in)=y
 



Functions for (multi)tapering
------------------------------------

multitap
-------
	[f, s] = multitap(x, taper)
	%%x a data vector, of length N; taper a matrix (N-by-p) of p tapers
	%%each of length N, as the data to be tapered
	%%returns the frequency f normalized to be in [0, 0.5],  and the
	%%multitapered spectrum estimate s, as the average of the p individual
	%%eigenspectra.  the individual eigenspectra are just the modulus
	%%square of the fft of the data x multiplied by one of the tapers.
	%%plots the spectrum on Decibel scale (10*log10(s)) vs. frequency
	%%no taper at all, the periodogram, corresponds to taper=ones(N,1);

	


EXAMPLE 
-------

Produce the periodogram and a 6-general prolate multitapered spectrum
estimate for the first N=2048 observations of GONG Month 10, l=85, m=0
data.  The values are in the file called "data", and the corresponding
gap structure in the file "Igaps".  Following are the matlab commands.


	load Igaps; I=Igaps; load data; N=2048; p=6; s=3; w=4/N;
	[x, m]=initprol(N, p, I, w);
	[prolval, prolvec]=prolate(x, m, I, p, N, s, w);
	prolvec=sqrt(N)*prolvec;
	taper=ones(N,1);
	[f, s]=multitap(data, taper);		%%periodogram
	[f1, s1]=multitap(data, prolvec);	%%6 gen.prol.tap. estim.
	plot(f, 10*log10(s),'r');
	hold
	plot(f1, 10*log10(s1),'g')
	title('Peridogarm (red) and 6 Multitapered Estimate (green)');
	xlabel('Frequency (cycles/unit time)');
	ylabel('Power (dB)');
	axis([0.1 0.3 -20 50]);			%%frequency band
	



*********************************************************	
*							*
*	CMEX FILES       				*
*							*
*		matlab functions that call		*
*		the C routines           		*
*********************************************************



GENERAL PROLATE 


initprol
--------
	Given the observation length N (integer), the number of tapers
	p (integer), the gap structure I (N-by-1 vector of zeroes and
	ones, I_i=1 if no gap at observation i, I_i=0 o.w.), and the
	frequency bandwidth w (real, measured in cycles/unit time, 
	0 < w < 0.5), initprol calculates the general sine tapers x
	(N-by-p vector, also used as initial guess vector for the 
	general prolate algorithm, see later), and the first row of
	the prolate Toeplitz matrix m (N-by-1 vector, to be used in the
	general prolate algorithm).  The initial guess vectors are
	orthonormal.

	Example:   N=32; p=2; w=0.125; I=ones(N,1);
		   [x, m]=initprol(N, p, I, w);


prolate
-------
	Given an initial guess vector x (N-by-p), the first row
	of the Toeplitz matrix m (N-by-1), the gap structure I 
	(N-by-1 vector os zeroes and ones, I_i=1 if no gap at 
	observation i, I_i=0 o.w.), the number of tapers p, 
	the length of the observations N, the number of Lanczos 
	vectors needed in the iterations s (s=3 suggested), and the
	bandwidth w (0 < w < 0.5), prolate estimates the p vectors,
	prolvec, with the largest energy concentration in the
	cycles/unit time frequency band [-w, w], along with the 
	concentrations, prolval.  The resulting tapers are
	orthonormal.  	

	Example:  N=32; p=2; w=0.125; I=ones(N,1);  s=3;
		  [x, m]=initprol(N, p, I, w);
		  [prolval, prolvec]=prolate(x, m, I, p, N, s, w);
		  prolvec=sqrt(N)*prolvec;





GENERAL MINIMUM BIAS

initmb
--------
	
	Given the observation length N (integer), the number of tapers
	p (integer), the gap structure I (N-by-1 vector of zeroes and
	ones, I_i=1 if no gap at observation i, I_i=0 o.w.), 
	initmb calculates the general sine tapers x
	(N-by-p vector, also used as initial guess vector for the 
	general minimum bias algorithm, see later), and the first row of
	the minimum bias Toeplitz matrix m (N-by-1 vector, to be used
	in the general minimum bias algorithm).  

	Example:   N=32; p=2; I=ones(N,1);
		   [x, m]=initmb(N, p, I);


minbias
-------
	Given an initial guess vector x (N-by-p), the first row
	of the Toeplitz matrix m (N-by-1), the gap structure I 
	(N-by-1 vector os zeroes and ones, I_i=1 if no gap at 
	observation i, I_i=0 o.w.), the number of tapers p, 
	the length of the observations N, the number of Lanczos 
	vectors needed in the iterations s (s=3 suggested),
	minbias estimates the p least locally biased vectors,
	mbvec, and their associated local biases (1-mbval).
	The resulting tapers are orthonormal.

	Example:  N=32; p=2; I=ones(N,1);  s=3;
		  [x, m]=initmb(N, p, I);
		  [mbval, mbvec]=minbias(x, m, I, p, N, s);
		  mbvec=sqrt(N)*mbvec;


	



*********************************************************	
*							*
*	C EXECUTABLES    				*
*							*
*********************************************************



GENERAL PROLATE

init.prol
---------
	
	Given the observation length N (integer), the number of tapers
	p (integer), the frequency bandwidth w (real, measured in
	cycles/unit time, 0 < w < 0.5), and the gap structure I
	(N-by-1 vector of zeroes and ones, I_i=1 if no gap at
	observation i, I_i=0 o.w.), init.prol calculates the general
	sine tapers x (N-by-p vector, also used as initial guess
	vector for the general prolate algorithm, see later), and the
	first row of the prolate Toeplitz matrix m (N-by-1 vector, to
	be used in the general prolate algorithm).  The initial guess
	vectors are orthonormal.

	Prompts the user to enter N, p, w, filename containing I, 
	filename to save the initial guess vector (general sine
	taper), and filename to save the first row of the general 
	prolate	Toeplitz matrix.  Filenames are character vectors of
	length not to exceed 20 characters.  The gap structure, I, is
	supposed to be an ASCII file; the results are saved into 
	binary files, say file1 and file2, as doubles.  To read them
	into matlab, and then properly scale them:
		
		% matlab
		>> fid=fopen('file1','r');
		>> x=fread(fid, [N,p], 'double');
		>> x=sqrt(N)*x;
		
prol
----
	Given an initial guess vector x (N-by-p), the first row
	of the Toeplitz matrix m (N-by-1), the gap structure I 
	(N-by-1 vector of zeroes and ones, I_i=1 if no gap at 
	observation i, I_i=0 o.w.), the number of tapers p, 
	the length of the observations N, the number of Lanczos 
	vectors needed in the iterations s (s=3 suggested), and the
	bandwidth w (0 < w < 0.5), prolate estimates the p vectors,
	prolvec, with the largest energy concentration in the
	cycles/unit time frequency band [-w, w], along with the 
	concentrations, prolval.  The resulting tapers are
	orthonormal. 
	
	Prompts for N, p, w, s, filename with I, filename with initial
	vector estimate, filename with first row of the Toeplitz
	matrix, filename to save prolval and prolvec.  All filenames
	must have less than 20 characters.  

	
	Example:  % init.prol
		  Enter N
		  32
		  Enter p
		  2
		  Enter the bandwidth
		  0.125
		  Enter the filename containing I
		  I
		  Enter the filename to save the initial guess vector
		  file1
		  Enter the filename to save row1 of Toeplitz
		  file2

	          % prol
		  Enter N
		  32
		  Enter p
		  2
		  Enter the bandwidth
		  0.125
		  enter the number of lanczos vectors needed
		  3
		  enter the filename containing I
		  I
		  enter the filename with initial vector estimates
		  file1
		  enter the filename containing row 1 of M
		  file2
		  enter the filename to save the estimated e.values
		  file3
		  enter the filename to save the estimated e.vectors
		  file4

	Read the vectors into matlab:

		% matlab
		>> fid=fopen('file3','r')
		>> conc=fread(fid, inf, 'double');	
		>> conc =
		   0.99998366170567		%%power concentration 
		   0.99965976582312		%%in band (-w, w).

		>> fid=fopen('file4','r')
		>> tap=fread(fid, [32, 2], 'double');
		>> tap=sqrt(32)*tap;
		>> plot(tap(:,1));		%%optimal taper
		>> plot(tap(:,2));		%%second best taper
prolfft
--------
	As prol, except that it uses an fft approximation to some
	matrix products, and thus saves time at slight precision 
	expense. 

		  % prolfft
		  Enter N
		  32
		  Enter p
		  2
		  Enter the bandwidth
		  0.125
		  enter the number of lanczos vectors needed
		  3
		  enter the filename containing I
		  I
		  enter the filename with initial vector estimates
		  file1
		  enter the filename containing row 1 of M
		  file2
		  enter the filename to save the estimated e.values
		  file5
		  enter the filename to save the estimated e.vectors
		  file6

		


GENERAL MINIMUM BIAS


init.mb		
-------
	Given the observation length N (integer), the number of tapers
	p (integer), the gap structure I (N-by-1 vector of zeroes and
	ones, I_i=1 if no gap at observation i, I_i=0 o.w.), 
	init.mb calculates the general sine tapers x (N-by-p vector,
	also used as initial guess vector for the general minimum bias
	algorithm, see later), and the first row of the minimum bias
	Toeplitz matrix m (N-by-1 vector, to be used in the general
	minimum bias algorithm).   

	Initial guess (general sine tapers) and first row of the 
	general minimum bias taper Toeplitz matrix.  It works as
	init.prol, except there is no need for the
	bandwidth parameter w.  
	
mb
---	Given N, p, s, I, initguess, and row1, mb estimates the first p
	general minimum bias tapers, and saves them into files
	specified by the user.  


mbfft
-----
	FFT approximation used in algortihm; faster but less precise
	than mb.


	Example:  % init.mb
		  Enter N
		  32
		  Enter p
		  2
		  Enter the filename containing I
		  I
		  Enter the filename to save the initial guess vector
		  file1
		  Enter the filename to save row1 of Toeplitz
		  file2

	          % mb
		  Enter N
		  32
		  Enter p
		  2
		  enter the number of lanczos vectors needed
		  3
		  enter the filename containing I
		  I
		  enter the filename with initial vector estimates
		  file1
		  enter the filename containing row 1 of M
		  file2
		  enter the filename to save the estimated e.values
		  file3
		  enter the filename to save the estimated e.vectors
		  file4

		  % mbfft
		  Enter N
		  32
		  Enter p
		  2
		  enter the number of lanczos vectors needed
		  3
		  enter the filename containing I
		  I
		  enter the filename with initial vector estimates
		  file1
		  enter the filename containing row 1 of M
		  file2
		  enter the filename to save the estimated e.values
		  file5
		  enter the filename to save the estimated e.vectors
		  file6

		  % matlab
		  >> fid=fopen('file4','r');
		  >> tap=fread(fid, [32,2],'double');
		  >> fid=fopen('file6','r');
		  >> tapfft=fread(fid, [32,2],'double');
		  >> plot(tap,'b'); hold
		  >> plot(tapfft,'r'); 
		  >> title('Tapers in blue, FFT approximated tapers in red')

		
	
EXAMPLE 
-------

Produce the periodogram and a 4-general prolate multitapered spectrum
estimate for the first N=2048 observations of GONG Month 10, l=85, m=0
data.  The values are in the ASCII file called "data", and the corresponding
gap structure in the file "Igaps".  Following are the matlab commands.

Example:  % init.prol
		  Enter N
		  2048
		  Enter p
		  4
		  Enter the bandwidth
		  0.002
		  Enter the filename containing I
		  Igaps
		  Enter the filename to save the initial guess vector
		  file1
		  Enter the filename to save row1 of Toeplitz
		  file2

	          % prolfft
		  Enter N
		  2048
		  Enter p
		  4
		  Enter the bandwidth
		  0.002
		  enter the number of lanczos vectors needed
		  3
		  enter the filename containing I

		  Igaps
		  enter the filename with initial vector estimates
		  file1
		  enter the filename containing row 1 of M
		  file2
		  enter the filename to save the estimated e.values
		  file3
		  enter the filename to save the estimated e.vectors
		  file4

		  % matlab
			>> load data
			>> N=2048; p=4
			>> fid=fopen('file4','r')
			>> prolvec=fread(fid, [N,p],'double');
			>> prolvec=sqrt(N)*prolvec;
			>> taper=ones(N,1);
			>> [f, s]=multitap(data, taper);	%%periodogram
			>> [f1, s1]=multitap(data, prolvec);	
						%%4 gen.prol.tap. estim.
			>> plot(f, 10*log10(s),'r');
			>> hold
			>> plot(f1, 10*log10(s1),'b')
			>> hold off
			>> title('Peridogarm (red) and 4 Multitapered
					Estimate (green)');
			>> xlabel('Frequency (cycles/unit time)');
			>> ylabel('Power (dB)');
			>> axis([0.1 0.3 -10 50])  %%detail