http://www.k-wave.org/forum/topic/fft_x_temp-and-temp_1_rs3d Support for the k-Wave MATLAB toolbox en-US Tue, 12 May 2026 23:29:44 +0000 http://bbpress.org/?v=1.0.2 q http://www.k-wave.org/forum/search.php http://www.k-wave.org/forum/topic/fft_x_temp-and-temp_1_rs3d#post-3735 Thu, 20 Jun 2013 22:23:13 +0000 Bradley Treeby 3735@http://www.k-wave.org/forum/ <p>Hi Anthony,</p> <p>The short answer is yes, but probably not for some time. We have done some preliminary work on this in 1D and it works well. However, the extension to 2D and 3D is still an active area of research.</p> <p>Brad. </p> http://www.k-wave.org/forum/topic/fft_x_temp-and-temp_1_rs3d#post-3729 Tue, 18 Jun 2013 16:39:31 +0000 Anthony 3729@http://www.k-wave.org/forum/ <p>Hi Brad,</p> <p>I just came across this topic and I am curious about the possibility to use non uniform grids: do you have any idea of when you plan to release such a feature (if it is still part of the plan)? It must be a tough job to do that!</p> <p>It would we awesome: I am simulating focused fields and I would love to be able to refine the grid near the focal point.</p> <p>Thanks a lot for what you and your team are doing with k-wave!<br /> Anthony</p> <p>PS : sorry if it not the right place to ask, but I didn't want to create a new topic only to satisfy my curiosity :-) </p> http://www.k-wave.org/forum/topic/fft_x_temp-and-temp_1_rs3d#post-1058 Fri, 23 Nov 2012 01:18:54 +0000 Bradley Treeby 1058@http://www.k-wave.org/forum/ <p>Hi Chao,</p> <p>Nonuniform in the sense that the grid points are not uniformly distributed, so dx/dy/dz are not constant across the domain. You can already define Nx/Ny/Nz and dx/dy/dx to not be equal if you wish, but the grid spacing is forced to be scalar.</p> <p>Brad. </p> http://www.k-wave.org/forum/topic/fft_x_temp-and-temp_1_rs3d#post-1057 Fri, 23 Nov 2012 00:17:10 +0000 huangchao 1057@http://www.k-wave.org/forum/ <p>Hi Brad</p> <p>Could you elaborate more on the nonuniform grid. Does nonuniform mean that the dx, dy and dz are not equal or Nx, Ny and Nz are not equal? Or it means dx (and/or dy, dz) is not a constant?</p> <p>Thanks,<br /> Chao </p> http://www.k-wave.org/forum/topic/fft_x_temp-and-temp_1_rs3d#post-1054 Thu, 22 Nov 2012 23:34:46 +0000 Bradley Treeby 1054@http://www.k-wave.org/forum/ <p>Hi Chao,</p> <p>The flag (and some additional functions) will be used in future releases to allow nonuniform grids to be defined. It's included in the current release so that HDF5 input files created with k-Wave version 1.0 are forward compatible. </p> <p>Brad. </p> http://www.k-wave.org/forum/topic/fft_x_temp-and-temp_1_rs3d#post-1052 Thu, 22 Nov 2012 21:59:11 +0000 huangchao 1052@http://www.k-wave.org/forum/ <p>Hi Jiri,</p> <p>Another quick question:</p> <p>I find in Appendix B of the manual that the 'nonuniform_grid_flag' must be set to 0. If so, what is the 'nonuniform_grid_flag' used for?</p> <p>Thanks,<br /> Chao </p> http://www.k-wave.org/forum/topic/fft_x_temp-and-temp_1_rs3d#post-1051 Thu, 22 Nov 2012 18:20:31 +0000 huangchao 1051@http://www.k-wave.org/forum/ <p>Hi Jiri,</p> <p>Thanks for your details explanation, I was confused by this statement:</p> <p>Get_FFT_X_temp().Compute_iFFT_3D_C2R(Get_Temp_1_RS3D());</p> <p>The reason is that I thought Get_Temp_1_RS3D() is the input argument, so I didn't understand where Get_Temp_1_RS3D get updated and wonder why iFFT still operates on Get_Temp_1_RS3D.</p> <p>I just realize Get_Temp_1_RS3D() is the output matrix defined in the FFTWComplexMatrix.h file, and it's very clear to me now.</p> <p>I am writing it here in case someone has similar problems as me.</p> <p>Thanks again!</p> <p>Chao </p> http://www.k-wave.org/forum/topic/fft_x_temp-and-temp_1_rs3d#post-1049 Thu, 22 Nov 2012 05:50:49 +0000 Jiri Jaros 1049@http://www.k-wave.org/forum/ <p>Hi Chao,</p> <p>Matrices <code>FFTW_X_temp, FFTW_Y_temp, and FFTW_Z_temp</code> serve as temporary storage for data in the complex domain at various places of the code. </p> <p>Any time I want to take FFT of something, I need to store the result somewhere. I usually need to take FFTs of three different spatial components, e.g. ux, uy, uz calculate duxdx, dyuduy, duzdz, so that’s why I need three temporary matrices. In the special case of pressure, I don’t need to take FFT of pressure 3 times, as the result will always be the same. That’s why I do so only once and store the result in <code>FFT_X_temp</code>. That’s the first line of <code>Compute_ddx_kappa_fft_p()</code>.</p> <p>However, when I multiply the FFT(p) with ddx/ddy/ddz and kappa I get 3 different results. Thus I need two more temporary matrices and store the output in FFT_X_temp, FFT_Y_temp and FFT_Z_temp in order to take the inverse FFT to calculate ux, uy, uz.</p> <p>Let’s dismantle the Matlab code for this part (the C++ code follows the MATLAB):</p> <pre><code>% calculate ux, uy and uz at the next time step using dp/dx, dp/dy and % dp/dz at the current time step ux_sgx = bsxfun(@times, pml_x_sgx, ... bsxfun(@times, pml_x_sgx, ux_sgx) ... - dt./rho0_sgx .* real(ifftn( ... bsxfun(@times, ddx_k_shift_pos, kappa .* p_k) ))); uy_sgy = bsxfun(@times, pml_y_sgy, ... bsxfun(@times, pml_y_sgy, uy_sgy) ... - dt./rho0_sgy .* real(ifftn( ... bsxfun(@times, ddy_k_shift_pos, kappa .* p_k) ))); uz_sgz = bsxfun(@times, pml_z_sgz, ... bsxfun(@times, pml_z_sgz, uz_sgz) ... - dt./rho0_sgz .* real(ifftn( ... bsxfun(@times, ddz_k_shift_pos, kappa .* p_k) )));</code></pre> <p>In order to save computational resources I calculate all of these statements simultaneously (kappa and p_k are the same think in all 3 statements).</p> <p>1) We take fft(p) and store the result in <code>p_k</code>. (in the C code this is<code>FFT_X_temp</code>)<br /> 2) We calculate</p> <p> a. <code>bsxfun(@times, ddx_k_shift_pos, kappa .* p_k)</code><br /> b. <code>bsxfun(@times, ddy_k_shift_pos, kappa .* p_k)</code><br /> c. <code>bsxfun(@times, ddz_k_shift_pos, kappa .* p_k)</code></p> <p>We’ll do this at once meaning taking an element from p_k, multiplying it with kappa (which is same in all statements) and store this in a single float variable say <code>p_k_element</code>. Then we carry out bsxfun on this element and ddx/ddy/ddz and store three results into target elements in <code>FFT_X_temp, FFT_Y_temp, FFT_Z_temp</code></p> <p>3) Now, we will take inverse FFTs and we need some place where to store it. We can’t override <code>ux_sgx</code> or <code>uy_sgy</code> or <code>uz_sgz</code> because we will need them. So we have to do</p> <p> a. <code>Temp_1_RS3D = real(ifft(FFT_X_temp));</code><br /> b. <code>Temp_2_RS3D = real(ifft(FFT_Y_temp));</code><br /> c. <code>Temp_3_RS3D = real(ifft(FFT_Z_temp));</code></p> <p>4) Now we finish the rest of the statements:</p> <p> a. <code>ux_sgx = bsxfun(@times, pml_x_sgx, bsxfun(@times, pml_x_sgx, ux_sgx) .- dt./rho0_sgx .* Temp_1_RS3D );</code><br /> b. <code>uy_sgy = bsxfun(@times, pml_y_sgy, bsxfun(@times, pml_y_sgy, uy_sgy) .- dt./rho0_sgy .* Temp_2_RS3D );</code><br /> c. <code>uz_sgz = bsxfun(@times, pml_z_sgz, bsxfun(@times, pml_z_sgz, uz_sgz) .- dt./rho0_sgz .* Temp_3_RS3D );</code> </p> <p>Please note that, the <code>Temp_1_RS3D, Temp_2_RS3D and Temp_3_RS3D</code> matrices are used at different places in the code in different meanings. Simply said, always I need a memory buffer I use one of these matrices since dynamic allocation could be very very very expensive!</p> <p>The code makes a lot of compromises to the readability in order to squeeze every possible FLOP and B/s from the HW :-)</p> <p>Jiri </p> http://www.k-wave.org/forum/topic/fft_x_temp-and-temp_1_rs3d#post-1048 Thu, 22 Nov 2012 04:53:04 +0000 huangchao 1048@http://www.k-wave.org/forum/ <p>Hi Jiri,</p> <p>I am studying the code and get confused by the temporary variables 'FFT_X_temp' and 'Temp_1_RS3D' (also for the temporary variables in the y and z directions).</p> <p>In the function 'Compute_ddx_kappa_fft_p', it seems that 'FFT_X_temp' is updated and used to stored the value of 'bsxfun(\@times, ddx_k_shift_pos, kappa .* p_k)', while in the function 'Compute_uxyz', the value of 'bsxfun(\@times, ddx_k_shift_pos, kappa .* p_k)' looks like being stored in 'Temp_1_RS3D', and I don't know where the 'Temp_1_RS3D' gets updated.</p> <p>Thanks,<br /> Chao </p>