% schts.m  --- 4 Sines/4 Cosines/Hartley Transforms
help schts
  S1=[];S2=[];S3=[];S4=[]; n=8
  C1=[];C2=[];C3=[];C4=[]; H=[]; 
for j=1:n
 for k=1:n
  S1(j,k) = sqrt(2/(n+1))*sin(pi*j*k/(n+1));   %-------------------Sine I
  S2(j,k) = sqrt(1/(2-mod(k,n)/k))*sqrt(2/n)*sin(pi*(j-1/2)*k/n); %Sine II
  S3(j,k) = sqrt(1/(2-mod(j,n)/j))*sqrt(2/n)*sin(pi*j*(k-1/2)/n); %Sine III
  S4(j,k) = sqrt(2/n)*sin(pi*(j-1/2)*(k-1/2)/n);                  %Sine IV

 C1(j,k) = sqrt((1+mod(1,j))/(n-1))*sqrt(1/(2-mod(1,k)))* ... % Cosine I
           sqrt(1/(2-mod(j,n)/j))*sqrt(1/(2-mod(k,n)/k))* ...
           cos(pi*(j-1)*(k-1)/ (n-1)); 
 C2(j,k) = sqrt((1+mod(1,k))/n)*cos(pi*(j-1/2)*(k-1)/n);  % Cosine II 
 C3(j,k) = sqrt((1+mod(1,j))/n)*cos(pi*(j-1)*(k-1/2)/n);  % III / NB C3=C2' 
 C4(j,k) = sqrt(2/n)*cos(pi*(j-1/2)*(k-1/2)/n);           % Cosine IV

  H(j,k) = sqrt(1/(n+1))*cos(2*pi*(j-1)*(k-1)/(n+1)) + ...  % Hartley
           sqrt(1/(n+1))*sin(2*pi*(j-1)*(k-1)/(n+1));
 end
end
  a=ones(n,1);A=full(spdiags([-a,2*a,-a], -1:1,n,n)); clear a
  A(1,1)=1; A(n,n)=1;    %---------------Test matrix for Cosine
  B=A; B=B-diag(diag(B));B=B+2*eye(n);%----Test matrix for Sine

 who S* C* H* A
 disp('Verify that (1)  C2''*C2 = eye(n) or Cj''*Cj or Sj''*Sj etc')
 disp('            (2)  C2''*A*C2 = diag ')
 disp('            (3)  S1''*B*S1 = diag ')
 
% ---- Chapter 4 
%  S2(j,k) = sqrt((1+mod(1,k))/n)*sin(pi*(j-1/2)*k/n); %Sine II - PM03 wrong

