%
%
%
% Optical Imaging and Spectroscopy
%
% David J. Brady
% Duke University
% www.opticalimaging.org
%
% Figure 3.17
%
%  spline representation
%
%
% First find spline function
%
clear all;
close all;
figure(1);set(gcf,'color','white');
for dip=1:2
    for m=0:3;
        xA=linspace(-7.485, 8.5, 1024);
        xS=xA(5)-xA(4);
        if m==0
            betaM=rect(xA-0.5);
        else
            betaM=conv(rect(xA+0.5),rect(xA-0.5));
            betaM=xS*betaM(481:1504);
            for cc=2:m
                betaM=conv(rect(xA-((-1)^cc)*0.5),betaM);
                betaM=xS*betaM(481:1504);
            end
        end
        %
        %
        % Find interpolating function
        %
        nC=20;
        p=zeros(1,nC);
        F=@(x)QQ(x,m);
        p(1)=quad(F,0,1);
        for k=2:nC
            F=@(x)(QQ(x,m).*cos(2*pi*(k-1)*x));
            p(k)=quad(F,0,1,10^-7);
        end
        checkP=p(1)*betaM;
        for shift=2:16
            checkP=checkP+p(shift)*(circshift(betaM',(shift-1)*64)'+circshift(betaM',-(shift-1)*64)');
        end
        figure(1);set(gcf,'color','white');
        if dip==2
            f=sin(pi*xA)+2.5+.8*cos(.3*pi*xA)+.001*xA.^3+5*exp(-20*(xA-pi).^2);
        else
            f=xA.^2/10;
        end
        %
        %
        % sampling matrix
        %
        ds=zeros(16,1024);
        us=zeros(16,1024);
        for pip=1:16
            ds(pip,:)=circshift(betaM',(8-pip)*64)';
            us(pip,:)=circshift(checkP',(8-pip)*64)';
        end
        fd=ds*f'/64;
        fu=fd'*ds;
        ms=num2str(m);
        subplot(4,2,2*m+dip);plot(xA,[f; fu],'-k');title(['m=' ms]);
        axis([-7.5 8.5 0 7]);
    end
end