Hi jlc4kk,
If you look at Fig. 3 in this paper, you can see that the velocity components each live on a staggered grid (ux is staggered in the x-direction, and so on). Equation 46 in the same paper then shows how to do the splitting, where using Einstein summation notation:
u_i^p = k_i k_j u_j
u_i^s = (delta_ij - k_i k_j) u_j
Because the split field components in each direction also depend on the velocity in other directions (which are staggered in different directions), these components must first be shifted before they are combined. The shift operators are labelled verbosely, so it should be clear what they do (shift_neg_y_pos_x will shift in the negative y-direction and positive x-direction).
In 2D, there are only two possibilities, moving from the x-staggered points to the y-staggered points, and vice versa. In 3D, there will be six possibilities (x <-> y, x <-> z, y <-> z), so you should create the required variables. Here, exp(1i*kgrid.kx*kgrid.dx/2) will produce a positive shift in the x-direction, while exp(-1i*kgrid.kx*kgrid.dx/2) will produce a negative shift (similarly for other directions).
It should then be clear how to calculate the various components. For example, the compression component of the particle velocity in the x-direction is given by
u_x^p = k_x k_x u_x + k_x k_y u_y + k_x k_z u_z
In matlab code, this would look like
ux_sgx^p = real(ifftn( kx_norm.^2 .* fftn(ux_sgx) + kx_norm .* ky_norm .* shift_neg_y_pos_x .* fftn(uy_sgy) + kx_norm .* kz_norm .* shift_neg_z_pos_x .* fftn(uz_sgz) ));
Hope that makes sense.
Brad.