function [psi] = psi_compute(T,q);
% Compute PSI Matrix for Cosine Transformation
% Inputs"
% T = time series length
% q = number of cosine transforms
% Output:
% psi = Txq matrix

  fvec = pi*(1:1:q);
  tvec = (0.5:1:T)'/T;
  psi = (sqrt(2)/T)*cos(tvec*fvec);
  fvec_2T = fvec/(2*T);
  ljt = (sin(fvec_2T))./fvec_2T;
  psi = psi.*repmat(ljt,T,1);
  
end

