%%
SLsample = false; % true; % false;

xlisthi = [0:127]*2*pi/128;

if SLsample
    im = phantom(128);
    flisthi = im(64,:);
else
    % define a function and its second derivative
    BLF = @(x) cos(3*x) + sin(5*x) - 2 + 0.35*cos(40*x);
    BLFdd = @(x) (-9*cos(3*x) - 25*sin(5*x)) -  0.35*40^2*cos(40*x)
    
    flisthi = BLF(xlisthi);
    fddlisthi = BLFdd(xlisthi);
end
figure(1); clf; hold on;

plot(xlisthi,50*flisthi,'k');
if SLsample
else
    plot(xlisthi,fddlisthi,'--k');
end
%% ----------  differentiate using convolution filter -----------------
f1 = [1 -2 1]   *128/(2*pi);
f2 = ([-1 16 -30  16 -1]/12) *128/(2*pi);
dxhi1 = (128/(2*pi))*conv(flisthi,f1,'same');
%dxhi1(1) = 0; % this prevents an edge effect.
%dxhi1 = dxhi1(2:129);
dxhi2 = (128/(2*pi))*conv(flisthi,f2,'same');
%dxhi2 = dxhi2(3:130);

% xz = (13+(0 - dxhi(13))/(dxhi(14) - dxhi(13)))*2*pi/128;
% 
% fxz = cos(3*xz) + sin(5*xz) -1;
% 

figure(1); plot(xlisthi,dxhi1,'r');

plot(xlisthi,dxhi2,'m+');

if SLsample
else
figure(3); clf ; hold on;
plot(xlisthi,fddlisthi - dxhi1,'r');

plot(xlisthi,fddlisthi - dxhi2,'m+');
end

%% ---------- differentiate in Fourier Domain ---------------

klist = [-64:63];

fh = fft(flisthi);
fhs = fftshift(fh);
figure(2); clf ; hold on;
figure(2);plot(klist,abs(fhs),'k');

figure(2);plot(klist,klist.^2/128,'--b'); 
dfs = fhs .* (-klist.^2);
plot(klist,abs(dfs./128),'r');
df = ifftshift(dfs);
d2 = ifft(df); % *2*pi/128;
 
figure(1); plot(xlisthi,d2,'--b');
if SLsample
else
figure(3);  plot(xlisthi,fddlisthi -d2,'--b');
end

%% --- show the second order derivative in spatial domain

sd2 = ifft( fftshift((-klist.^2)));
sd2 = ifftshift(sd2);

figure(4); plot(sd2);