% function TEMPORALBASESGEN constructs temporal basis vectors that spans
% the user-specified bandwidth of interest
%
% inputs: 1,2) [f1,f2] is the bandwidth of interest, e.g. for 1-30 Hz use f1 =1,
%         f2 = 30.
%         3) fs is the sampling frequency
%         4) T is the number of data time samples (as in [NxT] spatiotemporal
%         data)
%
% output: A [TxL] matrix C, where the L columns of C are the temporal basis vectors 
%
% Author: Tulaya Limpiti
% last updated: 2/23/2006

function C = TemporalBasesGen(f1,f2,fs,T)

% default error 
err = 0.99;

% covariance for bandpass white noise  
rw = 2*(f2-f1)/fs;
for M = 1:(T-1)
    rw = [rw, (sin(2*pi*M*f2/fs)-sin(2*pi*M*f1/fs))/(M*pi)];
end

Rc = toeplitz(rw);            
rank(Rc);

% eigendecomposition
L0 = 2*(f2-f1)*T/fs;
[V,D] = eig(Rc);
[D,iD] = sort(diag(D));
D = flipud(D);
iD = flipud(iD);
V = V(:,iD);                % sort eigenvector according to evals

% L = # eigenvalues such that er^2/tr(Rc) < (1-err)%, default = 1%
% er := tr(Rc) - tr(V^H*Rc*V);

L = 1; 
er = sum(D(1:L))/sum(D);

while(er <= err)
    L = L + 1;
    er = sum(D(1:L))/sum(D);
end

C = V(:,1:L);       % columns of C are temporal basis vectors