%% sample function
BLF = @(x) cos(5*x) + sin(3*x) - 1;

xlisthi = [0:127]*2*pi/128;
flisthi = BLF(xlisthi);% cos(5*xlisthi)  + sin(3*xlisthi)  - 1;

xlist1 = [0:3]*2*pi/4;
flist1 =  BLF(xlist1); % cos(5*xlist1)   + sin(3*xlist1) - 1;

xlist2 = [0:7]*2*pi/8;
flist2 = BLF(xlist2); % cos(5*xlist2)    + sin(3*xlist2)  - 1;

xlistnq = [0:9]*2*pi/10;
flistnq = BLF(xlistnq); % cos(5*xlistnq)  + sin(3*xlistnq) - 1;

xlist3 = [0:15]*2*pi/16;
%xlist3 = [0:12]*2*pi/12;
flist3 = BLF(xlist3); % cos(5*xlist3)+ sin(3*xlist3) - 1;

figure(1); clf; hold on;

plot(xlisthi,flisthi,'k');
plot(xlist2,flist2,'ob');
plot(xlist1,flist1,'+r');
plot(xlistnq,flistnq,'xm');
plot(xlist3,flist3,'*g');

%% ----------------------- Sinc Interpolation ---------------------
np =3;% number of repeated periods MUST BE ODD NUMBER.
per = 2*pi*[-floor(np/2):floor(np/2)];
xx1 =  reshape( (ones(np,1)*xlist1  + per' * ones(size(xlist1)))',1,[]);
xx2 =  reshape( (ones(np,1)*xlist2  + per' * ones(size(xlist2)))',1,[]);
xxnq =  reshape( (ones(np,1)*xlistnq  + per' * ones(size(xlistnq)))',1,[]);
xx3 =  reshape( (ones(np,1)*xlist3  + per' * ones(size(xlist3)))',1,[]);
xxhi = reshape( (ones(np,1)*xlisthi + per' * ones(size(xlisthi)))',1,[]);


nn = length(xxhi);
if1 = zeros(nn,1);
if2 = if1;
if3 = if1;
ifnq = if1;

for n = 1:nn
    tt1 = 2*(xx1 - xxhi(n));
    ftn1 = sinc(tt1/pi); % note Matlab Sinc scales by pi
 
    if1(n) = sum( reshape(flist1'*ones(1,np),[],1).* ftn1');
 
    tt2 = 4*(xx2 - xxhi(n));
    ftn2 = sinc(tt2/pi);

    if2(n) = sum( reshape(flist2'*ones(1,np),[],1).* ftn2');

     ttnq = 5*(xxnq - xxhi(n));
    ftnnq = sinc(ttnq/pi);

    ifnq(n) = sum( reshape(flistnq'*ones(1,np),[],1).* ftnnq');
   
    
    tt3 = 8*(xx3 - xxhi(n));    
%    tt3 = 6*(xx3 - xxhi(n));    
    ftn3 = sinc(tt3/pi);
 
    if3(n) = sum( reshape(flist3'*ones(1,np),[],1).* ftn3');
end
%%
figure(3); clf; hold on;

plot(xlisthi,flisthi,'k');
plot(xlist2,flist2,'ob');
plot(xlist1,flist1,'+r');
plot(xlistnq,flistnq,'xm');
plot(xlist3,flist3,'*g');

si = floor(np/2)*128;
ei = si+128;
plot(xlisthi,if1(si+1:ei),'--r');
plot(xlisthi,if2(si+1:ei),'--b');
plot(xlisthi,ifnq(si+1:ei),'--m');
plot(xlisthi,if3(si+1:ei),'g.');
title(['Sinc interpolation using ',num2str(np),' periods']);
%% ----------------------- inverse zero fill ----------------------
fx1 = fft(flist1)*128/4;
zfx1 = [fx1(1:2),zeros(1,124),fx1(3:4)];
izf1 = ifft(zfx1);

figure(1); plot(xlisthi,real(izf1),'--r')

fx2 = fft(flist2)*128/8;
zfx2 = [fx2(1:4),zeros(1,120),fx2(5:8)];
izf2 = ifft(zfx2);
 
figure(1); plot(xlisthi,real(izf2),'--b')


fxnq = fft(flistnq)*128/10;
zfxnq = [fxnq(1:5),zeros(1,118),fxnq(6:10)];
izfnq = ifft(zfxnq);
 
figure(1); plot(xlisthi,real(izfnq),'--m')

 
fx3 = fft(flist3)*128/16;
zfx3 = [fx3(1:8),zeros(1,112),fx3(9:16)];
%fx3 = fft(flist3)*128/12;
%zfx3 = [fx3(1:5),zeros(1,116),fx3(6:12)];
izf3 = ifft(zfx3);
 
 figure(1); plot(xlisthi,real(izf3),'g.')

title('Fourier zero fill interpolation');

%% -------------------plot Fourier ----------------------------------
figure(4);clf ; hold on;
ffh = fft(flisthi);
subplot(2,2,1); plot(xlisthi-pi,abs(fftshift(ffh)),'k');
subplot(2,2,2);plot(xlisthi-pi,abs(fftshift(zfx1)),'r');
subplot(2,2,3);plot(xlisthi-pi,abs(fftshift(zfx2)),'b');
subplot(2,2,4);plot(xlisthi-pi,abs(fftshift(zfxnq)),'m');
%subplot(2,2,4);plot(xlisthi-pi,abs(fftshift(zfx3)),'g');