% Compute Hermite Functions 

clear all
d=13;   %Odd numbers work best if N is divisble by 4 there is 1 repeated eigenvalue
        % and the eigenspace that Matlab returns for that value can have complex eigenvectors 
N=d^2; % Legacy for computing Gabor systems with these 

D = diag(ones(1,N),0) + diag(-ones(1,N-1),1);
D(N,1)=-1;      % D is circulant Diffrence Matrix;
T=eye(N)-D;
Id=eye(d);
L=(D'*D);       % Discrete laplacian
F=fft(eye(N))/sqrt(N);%DFT
X=(F*L*F');     %analog of potential function x^2 in Scrodinger operator- is diagonal
X=real(X);      %You can prove the diagonal of X is real see Maser Lammers paper for example
                % but for matlab combutations we take away any lingering
                % imaginary part for the diagonalization

[H,S,V]=svd(X+L);   % X+L is Gauss-Hermite Diff eq or Schrodinger operator
                    % use svd instead of eig to order the eignevectors 
                    % the eigenvalues are square roots of S
                    % F commutes with X+L so they share eigenvectors
                    % H orthogonal matrix
EF=H'*F*H;  %Show H diagonalizes F
figure(1)
plot((EF))  % plot illustrates it is diagoonal with values 1,-1,i and -i
figure(2)
imagesc(H*H')
H=-H; %sometimes returns -1 times Hermite functions depending on d
figure(3)
clf
plot((H(:,N-2:N))) %plot first 3 Hermite
figure(4)
clf
plot(T^(floor(N/2))*H(:,N-2:N))% Shifted to center
figure(5)
clf
fgauss=H(:,N);      %finite gausian

plot(fgauss-F*fgauss)   %compute difffernce of gaussian and the FFT of guassian. Some 
                        %small complex part leeks for the FFT of fgauss but esetially zero
figure(4)