k-Wave User Forum » Topic: k vector for even number of grid points

k-Wave User Forum » Topic: k vector for even number of grid points http://www.k-wave.org/forum/topic/k-vector-for-even-number-of-grid-points Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 06:51:06 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php mat.simurda on "k vector for even number of grid points" http://www.k-wave.org/forum/topic/k-vector-for-even-number-of-grid-points#post-4959 Fri, 23 Jan 2015 11:35:29 +0000 mat.simurda 4959@http://www.k-wave.org/forum/ Hi Brad, thank you very much for your reply. Yes, your explanation helped me to understand the problem. Have a nice day and once again, thank you. Best regards Matej Bradley Treeby on "k vector for even number of grid points" http://www.k-wave.org/forum/topic/k-vector-for-even-number-of-grid-points#post-4947 Thu, 22 Jan 2015 18:56:47 +0000 Bradley Treeby 4947@http://www.k-wave.org/forum/ 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 <a href="http://www.k-wave.org/papers/2011-Treeby-JASA.pdf">this paper</a>. Computationally you can check this: <pre><code>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))');</code></pre> 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: <pre><code>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');</code></pre> Hope that helps, Brad. mat.simurda on "k vector for even number of grid points" http://www.k-wave.org/forum/topic/k-vector-for-even-number-of-grid-points#post-4933 Thu, 08 Jan 2015 15:04:14 +0000 mat.simurda 4933@http://www.k-wave.org/forum/ Good day, thank you very much for the impressive K-wave toolbox. I believe I have found an error in your source code: While defining wave vector kgrid.kx_vec/ky_vec/kz_vec for even number of grid points you don't set the wave number for index -N/2 to zero. ... if rem(Nx, 2) == 0 % grid dimension has an even number of points nx = ((-Nx/2:Nx/2-1)/Nx).'; else ... This treats the highest wavenumber asymmetrically (see L.N. Trefethen, Spectral methods in MATLAB, Chapter 3). Best regards Matej