clc; clear; close all;

n = 16; 
h = 2*pi/n; 

xs = h*(1:n)';

%---- spectral differenciation

c = cos(2*xs) + cos(xs); dc_true = -2*sin(2*xs) - sin(xs);
% c = sin(2*xs) + sin(xs); dc_true = 2*cos(2*xs) + cos(xs);
% c = exp(-30*xs.^2); dc_true = -60*xs.*c;
% c = max(0, 1-abs(xs)); dc_true = (xs>=-1) - 2*(xs>=0) + (xs>=1);

% differenciation in Fourier domain
dc1 = real(ifft(1i*[0:n/2-1 0 -n/2+1:-1]'.*fft(c)));
% differenciation by Toeplitz matrix
dSN = [0 -.5*(-1).^(1:n-1).*cot((1:n-1)*h/2)];
tdSN = toeplitz([dSN(1);dSN(end:-1:2)'], dSN);
dc2 = tdSN*c;
% differenciation convolution in Fourier domain
dc3 = real(ifft(fft(dSN)'.*fft(c)));
% two-point differenciation
diff2pt = [0 1 zeros(1,n-3) -1]/2/h;
tdiff2pt = toeplitz([diff2pt(1);diff2pt(end:-1:2)'], diff2pt);
dc4 = tdiff2pt*c;

figure;
plot(xs,c,'b', xs,dc1,'r', xs,dc2,'g+', xs,dc3,'ro', xs,dc_true,'k');

figure; 
plot(xs,c,'b', xs,dc1,'r', xs,dc4,'g',xs,dc_true,'k');

%---- hilbert transform

% c = cos(2*xs) + cos(xs); hc_true = sin(2*xs) + sin(xs);
c = sin(2*xs) + sin(xs); hc_true = -cos(2*xs) - cos(xs);

hc = real(ifft((1/1i)*[ones(1,n/2) 0 -1*ones(1,n/2-1)]'.*fft(c)));

figure;
plot(xs,c,'b', xs,hc,'r', xs,hc_true,'k');

%---- spectral shift
dx = pi/4;

% c = cos(xs); tc_true = cos(xs-dx);
% c = cos(2*xs) + cos(xs); tc_true = cos(2*(xs-dx)) + cos(xs-dx);
% c = sin(2*xs) + sin(xs); tc_true = sin(2*(xs-dx)) + sin(xs-dx);
% c = exp(-30*(xs-pi).^2); tc_true = exp(-30*(xs-pi-dx).^2);
c = max(0, 1-abs(xs-pi)); tc_true = max(0, 1-abs(xs-pi-dx));

txs = xs-dx;

col = (2*pi)/h * sin( pi*txs/h ) ./ tan( txs/2);
t = toeplitz([col(1);col(end:-1:2)], col);

tc = t*c/n/n;

figure;
plot(xs,c,'b', xs,tc,'r', xs,tc_true,'k');

