Lagrange, idht, flt

Analysis/HermFour.m

@@ -1,5 +1,5 @@
 function [] = HermFour(f, n, m)
-[x, w]=GaussHermite(n); ww=log(sqrt(2)*w);
+[x, w]=GaussHermite(n,0,1); ww=log(sqrt(2)*w);
 H=zeros(m, n);
 H(1,:)=pi^(-1/4)*exp(-x.^2);
 H(2,:)=H(1,:).*(2*x);
@@ -21,6 +21,7 @@
     sum_log_scale=sum_log_scale+log_scale;
 end
 
+
 F=H'*diag(1i.^(0:m-1))*W;
 G=conj(F);
 [kx, ky]=meshgrid(sqrt(2)*x);

Analysis/sincFT2.m

@@ -10,9 +10,9 @@
 
 z=f(xx, yy);
 A=h^2/(2*pi);
-u=fftshift(fft2(fftshift(z)).*A); 
+u=fftshift(fft2(fftshift(z))*A); 
 %uu=-u./(kx.^2+ky.^2); uu(N/2+1,N/2+1)=0;
-v=fftshift(ifft2(fftshift(z))./A);
+v=fftshift(ifft2(fftshift(z))/A);
 
 
 figure(1); colormap(jet(256));

DPT/flt.m

@@ -0,0 +1,24 @@
+function [] = flt(n)
+% Fast Legendre Transform, based on Alpert-Rokhlin
+
+% Initialization
+k=32;
+s=4*k;
+
+% Step 1
+r=0:k-1;
+t=(1-cos((r+0.5)*pi/k))/2;
+tt=[t/2, (1+t)/2];
+
+% Step 2
+T=repmat(t,[k,1]);
+den=prod(T-T'+eye(k));
+l=0:s-1;
+U=repmat(l/s,[k,1])-repmat(t(:),[1,s]);
+U=bsxfun(@rdivide, bsxfun(@rdivide, prod(U), U), den(:));
+
+%Step 3
+h=log2(n/s)-1;
+
+end
+

DPT/fpt.m

@@ -1,5 +1,5 @@
 function [f] = fpt(a,b,c,g)
-% Fast Polynomial Transform, based on Potts (1998)
+% Fast Polynomial Transform, based on Potts (1998).
 % Basis exchange for orthogonal polynomials that satisfy the three term
 % recurrence relation: P[n](x)=(a[n]x+b[n])P[n-1](x)+c[n]P[n-2](x).
 % The last input argument, g, should be the coeffients of the linear
@@ -14,8 +14,8 @@
 
 global fileID;
 global create;
-
-N=length(g)-1; t=log2(N);
+N=length(g)-1;
+t=log2(N);
 
 filename=sprintf('fptU%d.bin',N);
 create=1-exist(filename,'file')/2;
@@ -25,7 +25,8 @@
     fileID=fopen(filename,'r');
 end
 
-gg=[g(:), zeros(N+1,1)]; gg(N+1,:)=[];
+gg=[g(:), zeros(N+1,1)];
+gg(N+1,:)=[];
 gg(N-1,1)=g(N-1)+c(N)*g(N+1);
 gg(N,:)=[g(N)+b(N)*g(N+1), a(N)*g(N+1)];
 ggg=zeros(N,2);
@@ -38,15 +39,15 @@
        j=2*i-1;
        h=m*(i-1)+1;
        [u11,u12,u21,u22]=calculateU(a,b,c,m-1,h,x);
-       ggg( i , : ) = c(h+1)*innerP(u11, u12, gg(j+2,:), gg(j+3,:));
+       ggg(i  , : ) = c(h+1)*innerP(u11, u12, gg(j+2,:), gg(j+3,:));
        ggg(i+1, : ) =        innerP(u21, u22, gg(j+2,:), gg(j+3,:));
-       ggg( i ,1:m) = ggg( i ,1:m)+gg(j,:);
+       ggg(i  ,1:m) = ggg(i  ,1:m)+gg(j  ,:);
        ggg(i+1,1:m) = ggg(i+1,1:m)+gg(j+1,:);
    end
    gg=ggg;
 end
 
-% f=gg(1,:)+(b(1)I+a(1)T')*gg(2,:)
+% f=gg(1,:)+(b(1)*I+a(1)*T')*gg(2,:)
 f=zeros(size(g)); 
 f(1)=gg(1,1)+b(1)*gg(2,1)+a(1)*gg(2,2)/2;
 f(2)=gg(1,2)+b(1)*gg(2,2)+a(1)*(gg(2,1)+gg(2,3)/2);
@@ -57,6 +58,7 @@
 end
 
 function [u11, u12, u21, u22] = calculateU(a,b,c,n,j,x)
+% Calculates the 4x4 matrix with polynomial entries
 global fileID;
 global create;
 if(create)
@@ -65,8 +67,7 @@
     fwrite(fileID, [u11; u12; u21; u22], 'double');
 else
     U=fread(fileID, [4, length(x)], 'double');
-    u11=U(1,:); u12=U(2,:);
-    u21=U(3,:); u22=U(4,:);
+    u11=U(1,:); u12=U(2,:); u21=U(3,:); u22=U(4,:);
 end
 end
 
@@ -84,6 +85,6 @@
 % Computes polynomial vector dot product [u1; u2]'*[a1; a2]
 a1(1)=a1(1)*sqrt(2); 
 a2(1)=a2(1)*sqrt(2);
-P=dct(u1.*idct(a1, length(u1))+u2.*idct(a2, length(u2)));
+P=dct(u1.*idct(a1,length(u1))+u2.*idct(a2,length(u2)));
 P(1)=P(1)/sqrt(2);
 end

DPT/idht.m

@@ -0,0 +1,11 @@
+function [z] = idht(f, N)
+% Inverse Discrete Hermite Transform
+x=sqrt(pi/(2*N))*(-N:2:N-2);
+H=zeros(N,N);
+H(1,:)=pi^-(1/4)*exp(-x.^2/2);
+H(2,:)=sqrt(2)*x.*H(1,:);
+for i=1:N-2
+    H(i+2,:)=sqrt(2/(i+1))*x.*H(i+1,:)-sqrt(i/(i+1))*H(i,:);
+end
+z=sqrt(2*pi/N)*H*f(x');
+end

Integration/GaussChebyshev.m

@@ -1,7 +1,7 @@
 function [x, w] = GaussChebyshev(a, b, n)
 % Returns abscissas and weights for the Gauss-Chebyshev n-point quadrature
-% over the interval [a, b]. Weigths include reciprocal weight function.
-th=((1:n)-1/2)*pi/n;
+% over the interval [a, b].
+th=pi*(1:2:2*n-1)/(2*n);
 x=(b-a)/2*cos(th)+(a+b)/2;
-w=(b-a)*pi/(2*n)*sin(th);
+w(1:n)=(b-a)*pi/(2*n);
 end

Integration/GaussHermite.m

@@ -1,7 +1,8 @@
-function [x, w] = GaussHermite(n)
+function [x, w] = GaussHermite(n, mu, sigma)
 % Returns abscissas and weights for the Gauss-Hermite n-point quadrature 
 % over the interval [-inf, inf] using the Golub-Welsch Algorithm.
 beta=sqrt((1:n-1)/2);
 [x,V]=trideigs(zeros(1,n), beta); 
-x=x'; w=sqrt(pi)*V(1,:).^2;
+x=(sqrt(2)*sigma)*x'+mu;
+w=sqrt(2*pi)*sigma*V(1,:).^2;
 end

Interpolation/Lagrange.m

@@ -1,20 +1,10 @@
 function [p] = Lagrange(x, y, t)
 % Evaluates the Lagrange interpolation polynomial.
-n=length(x);
-w=ones(size(x));
-for i=1:n
-    for j=1:n
-        if j~=i
-            w(i)=w(i)*(x(i)-x(j)); 
-        end  
-    end
-end
-g=1;
-p=0;
-for k=1:n
-    d=t-x(k);
-    p=p+y(k)./(w(k)*d);
-    g=g.*d;
-end
-p=g.*p;
+n=length(x); 
+m=length(t);
+X=repmat(x(:).',[n,1]);
+w=prod(X-X.'+eye(n));
+r=y./w;
+D=repmat(t(:).',[n,1])-repmat(x(:),[1,m]);
+p=prod(D).*(r*(1./D));
 end

Special Functions/HermitePsi.m

@@ -0,0 +1,12 @@
+function [y]=HermitePsi(a, x)
+% Evaluates the Hermite function series given by the coefficients a(n) 
+n=length(a);
+y=(n>1)*a(n); yy=zeros(size(x));
+for k=n-2:-1:1
+    temp=y;
+    y=a(k+1)+sqrt(2/(k+1))*x.*y-sqrt((k+1)/(k+2))*yy;
+    yy=temp;
+end
+h0=pi^(-1/4)*exp(-x.^2/2);
+y=h0.*(a(1)+sqrt(2)*x.*y-sqrt(1/2)*yy);
+end

Special Functions/jinc.m

@@ -0,0 +1,5 @@
+function [j] = jinc(a,x)
+t=(x==0);
+j=(besselj(a,x)+t)./(x+t);
+end
+

Special Functions/sinc.m

@@ -0,0 +1,5 @@
+function [y] = sinc(x)
+t=(x==0);
+y=(sin(x)+t)./(x+t);
+end
+