Hi Matej,
Thanks for your post. If you define the wavenumbers asymmetrically, the band limited interpolant (BLI) will have both real and complex parts (this is point being made in the caption of Fig. 3.1 in the book you mention). The exact form of the BLI is given in Eq. 28 in this paper.
Computationally you can check this:
clear all
% define x grid
dx = 1;
Nx = 20;
x = (0:Nx-1);
% define delta function
f = zeros(1, Nx);
f(1) = 1;
% define set of wavenumbers
k = (-pi/dx) : (2*pi/(dx*Nx)) : (pi/dx-2*pi/(dx*Nx));
% define fine grid for displaying BLI
x_fine = 0:dx/10:(Nx-1)*dx;
% calculate Fourier coefficients
func_k = ifftshift(fft(f)) / Nx;
% calculate BLI by summing sinusoids
BLI_fft = sum( func_k * exp(-1i * k.' * x_fine), 1);
% calculate BLI using closed form from JASA paper
BLI_exact = 1/Nx * (1i + cot(x_fine*pi/(dx*Nx))) .* sin(x_fine*pi/dx);
% plot BLI
figure;
subplot(2, 1, 1)
plot(x_fine, real(BLI_fft), 'k-');
hold on;
plot(x_fine, real(BLI_exact), 'b.');
plot(x, f, 'ro');
xlabel('x');
ylabel('real(f(x))');
title('Unsymmetric wavenumbers');
subplot(2, 1, 2)
plot(x_fine, imag(BLI_fft), 'k-');
hold on;
plot(x_fine, imag(BLI_exact), 'b.');
plot(x, zeros(size(x)), 'ro');
xlabel('x');
ylabel('imag(f(x))');
In this example, the BLI is correct in that it matches the function exactly at all of the grid points. However, there is a sinusoidal variation in the complex part that would give the derivative incorrectly. However, from Trefethen's example in Fig. 3.1, you'll note that taking the real part of the interpolant exp(i*pi*x/h) gives the correct result. This is exactly what happens in k-Wave - the real part is taken after each spectral derivative is computed.
Again, you can check this gives the same thing computationally:
clear all;
% define grid parameters (Nx must be even)
Nx = 10;
dx = 2*pi/Nx;
x = (0:Nx-1)*dx;
% define function
% 1: sinusoid
% 2: sawtooth
func = 1;
switch func
case 1
% sinusoid
f = sin(4*x);
dfdx = 4*cos(4*x);
case 2
% saw tooth
f = ones(1, Nx);
f(2:2:end) = 0;
dfdx = 0;
end
% define wavenumbers as they are in k-Wave
kx = (2*pi/dx) .* (-Nx/2:Nx/2-1)/Nx;
% put wavenumbers in correct order for FFT
kx = ifftshift(kx);
% compute spectral derivative
dfdx_spec = real(ifft(1i*kx .* fft(f)));
% define Trefethyn wavenumbers (this sets the Nyquist value to zero)
kx_tref = [0:Nx/2-1 0 -Nx/2+1:-1];
% compute spectral derivative
dfdx_tref = ifft(1i*kx_tref .* fft(f));
% plot derivative
figure;
subplot(2, 1, 1);
plot(x, f, 'k--');
hold on;
plot(x, dfdx, 'k-');
plot(x, dfdx_spec, 'b.');
plot(x, dfdx_tref, 'rd');
legend('f', 'dfdx - analytical', 'dfdx - k-Wave', 'dfdx - Trefethyn');
xlabel('x');
ylabel('f');
subplot(2, 1, 2);
plot(x, dfdx - dfdx_spec, 'b-');
hold on;
plot(x, dfdx - dfdx_tref, 'r--');
xlabel('x');
ylabel('error');
legend('k-Wave', 'Trefethyn');
Hope that helps,
Brad.