http://www.k-wave.org/forum/topic/timing-error-in-the-generated-rf-data Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 08:18:31 +0000 http://bbpress.org/?v=1.0.2 q http://www.k-wave.org/forum/search.php http://www.k-wave.org/forum/topic/timing-error-in-the-generated-rf-data#post-7558 Sat, 06 Jun 2020 14:36:47 +0000 Bradley Treeby 7558@http://www.k-wave.org/forum/ <p>Hi Wonseok,</p> <p>Could you try setting <code>medium.alpha_mode = &#39;no_dispersion&#39;</code>?</p> <p>When you have power law acoustic absorption, causality dictates that you must also have dispersion (a dependence of the sound speed on frequency). As your power law is very close to 1, I suspect that this is giving you a noticeably different sound speed in the medium. </p> <p>Have a look at the <a href="http://www.k-wave.org/documentation/example_na_modelling_absorption.php">Modelling Power Law Absorption Example</a> in MATLAB. You can see for the same background sound speed, the actual wave speed can vary quite a lot (in accordance with the theory) depending on the absorption parameters.</p> <p>Brad. </p> http://www.k-wave.org/forum/topic/timing-error-in-the-generated-rf-data#post-7522 Fri, 22 May 2020 03:29:36 +0000 Wonseok Choi 7522@http://www.k-wave.org/forum/ <p>Thank you for you reply.<br /> I have posted the working code.</p> <p>I tried removing the sensor.frequency_response line, but it didn't work either. </p> <pre><code>clearvars; %% ========================================================================= % SIMULATION % ========================================================================= % Environment setting PML_size = 40; % size of the PML in grid points --&gt; PML = perfectly matched layer. assumed to simulate open boundaries Nz = 1600; % number of grid points in the x (row) direction Nx = 1600; % number of grid points in the y (column) direction dz = 1e-4; % grid point spacing in the z direction [m] dx = 1e-4; % grid point spacing in the x direction [m] Xsize = Nx * dx; Zsize = Nz*dz; % define a binary line sensor sensor.mask = zeros(Nz, Nx); posXgrid = 613:3:994; posZgrid = 56; for idxEle = 1:length(posXgrid) sensor.mask(posZgrid, posXgrid(idxEle)) = 1; end % Temporal arguments SoundSpeed = 1510; SampleFreq = 50e6; dt = 1/SampleFreq; Nt = ceil(1.1 * sqrt(Xsize^2+Zsize^2)/SoundSpeed/dt)/2; % Transducer arguments transducer.pitch = 0.298e-3; % [m] transducer.numElement = 128; % [m] transducer.centerFreq = 5e6; % [1/sec] transducer.BW = 60; % [%] sensor.frequency_response = [transducer.centerFreq, transducer.BW]; % create the computational grid kgrid = kWaveGrid(Nz, dz, Nx, dx); kgrid.setTime(Nt, dt); % Medium Parameters medium.alpha_coeff = 0.75; % [dB/(MHz^y cm)] medium.alpha_power = 1.01; medium.sound_speed = SoundSpeed; % [m/s] % create initial pressure distribution using makeDisc NumSource_x = 3; NumSource_z = 3; Disc.magnitude = 1; % [Pa] Disc.radius = 2; % scaled by grid point spacing. numTarget = 9; if numTarget == 9 % 9 target disc = 0; for idx_x = 1:NumSource_x for idx_z = 1:NumSource_z Disc.z_pos = 384*(idx_z-2)/(NumSource_z+1)+Nz/2; % [grid points] Disc.x_pos = 384*(idx_x-2)/(NumSource_x+1)+Nx/2; % [grid points] temp = Disc.magnitude * makeDisc(Nz, Nx, Disc.z_pos, Disc.x_pos, Disc.radius); disc = disc + temp; end end elseif numTarget == 3 % 3 target diagonal disc = 0; for idx_x = 1:NumSource_x Disc.z_pos = 384*(idx_x-2)/(NumSource_z+1)+Nz/2; % [grid points] Disc.x_pos = 384*(idx_x-2)/(NumSource_x+1)+Nx/2; % [grid points] temp = Disc.magnitude * makeDisc(Nz, Nx, Disc.z_pos, Disc.x_pos, Disc.radius); disc = disc + temp; end end figure; imagesc(disc + sensor.mask); axis image; % smooth the initial pressure distribution and restore the magnitude source.p0 = smooth(kgrid, disc, true); % set the input arguements: force the PML to be outside the computational % grid; switch off p0 smoothing within kspaceFirstOrder2D input_args = {&#39;PMLInside&#39;, false, &#39;PMLSize&#39;, PML_size, &#39;PlotPML&#39;, false, &#39;Smooth&#39;, false, &#39;DataCast&#39;, &#39;gpuArray-single&#39;}; %% RF generation % run the simulation sensor_data_raw = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:}); %% Time reversal source2 = source; source2.p0=0; sensor2 = sensor; sensor2.time_reversal_boundary_data = sensor_data_raw; kgrid2 = kgrid; medium2 = medium; medium2.alpha_coeff = 0; p0_tr = kspaceFirstOrder2D(kgrid2, medium2, source2, sensor2, input_args{:}); figure; imagesc(p0_tr);</code></pre> http://www.k-wave.org/forum/topic/timing-error-in-the-generated-rf-data#post-7512 Thu, 21 May 2020 11:09:31 +0000 Bradley Treeby 7512@http://www.k-wave.org/forum/ <p>Hi Wonseok,</p> <p>Can you post a working code? I'll take a look. Does it work if you remove the <code>sensor.frequency_response = ...</code> line?</p> <p>Brad. </p> http://www.k-wave.org/forum/topic/timing-error-in-the-generated-rf-data#post-7506 Tue, 19 May 2020 06:56:58 +0000 Wonseok Choi 7506@http://www.k-wave.org/forum/ <p>Dear developers,<br /> Thank you so much for providing this wonderful simulator!</p> <p>We are trying to generate RF data for photoacoustic computed tomography that uses multiple linear arrays circularly positioned around the target.<br /> We tried the code below and performed time reversal for image reconstruction.<br /> However, we got a wrong reconstruction that looks like the time reversing was done a little bit more.<br /> (i.e. circular targets appearing as donut shapes)<br /> We then tried to reconstruct the image using self-developed delay-and-sum algorithm, and we still got a similar result with time reversal.</p> <p>We found an unexpected error in the RF data: the wavefronts appearing at earlier timings than expected.<br /> Using point source targets, we calculated the expected timing of the wavefront, and we found that the wavefronts in the RF data appeared at about 88% of the calculated delay.<br /> When we manually multiplied 0.88 at the calculated delays for delay-and-sum, the image was well reconstructed.<br /> We think time reversal operation was working all right, but the forward process was somehow not right. </p> <p>We have provided the code as the following lines.<br /> Sorry that we couldn't include the element positions here for you to test the code.<br /> If you need to execute the code, you can simply put a single 128-element linear array on the top of the grid.<br /> (We also tried that but it also gave the same error.) </p> <p>Thank you in advance for kindly taking your time on this issue.<br /> Wonseok Choi.</p> <pre><code>% Environment setting PML_size = 40; % size of the PML in grid points --&gt; PML = perfectly matched layer. assumed to simulate open boundaries Nz = 1600; % number of grid points in the x (row) direction Nx = 1600; % number of grid points in the y (column) direction dz = 1e-4; % grid point spacing in the z direction [m] dx = 1e-4; % grid point spacing in the x direction [m] Xsize = Nx * dx; Zsize = Nz*dz; % Detector position load(&#39;ElementPosition.mat&#39;) % define a binary line sensor sensor.mask = zeros(Nz, Nx); posX = pos_ele(:,1); posXgrid = round((posX-(max(posX)+min(posX))/2)/abs(dx)) + Nx/2; posZ = pos_ele(:,2); posZgrid = round((posZ-(max(posZ)+min(posZ))/2)/abs(dz)) + Nz/2; for idxEle = 1:length(posX)/8 sensor.mask(posZgrid(idxEle), posXgrid(idxEle)) = 1; end % Temporal arguments SoundSpeed = 1510; SampleFreq = 50e6; dt = 1/SampleFreq; Nt = ceil(1.1 * sqrt(Xsize^2+Zsize^2)/SoundSpeed/dt)/2; % Transducer arguments transducer.pitch = 0.298e-3; % [m] transducer.numElement = 128; % [m] transducer.centerFreq = 5e6; % [Hz] transducer.BW = 60; % [%] sensor.frequency_response = [transducer.centerFreq, transducer.BW]; % create the computational grid kgrid = kWaveGrid(Nz, dz, Nx, dx); kgrid.setTime(Nt, dt); % Medium Parameters medium.alpha_coeff = 0.75; % [dB/(MHz^y cm)] medium.alpha_power = 1.01; medium.sound_speed = SoundSpeed; % [m/s] % create initial pressure distribution using makeDisc NumSource_x = 3; NumSource_z = 3; Disc.magnitude = 1; % [Pa] Disc.radius = 2; % scaled by grid point spacing. numTarget = 9; disc = 0; for idx_x = 1:NumSource_x for idx_z = 1:NumSource_z Disc.z_pos = 384*(idx_z-2)/(NumSource_z+1)+Nz/2; % [grid points] Disc.x_pos = 384*(idx_x-2)/(NumSource_x+1)+Nx/2; % [grid points] temp = Disc.magnitude * makeDisc(Nz, Nx, Disc.z_pos, Disc.x_pos, Disc.radius); disc = disc + temp; end end % smooth the initial pressure distribution and restore the magnitude source.p0 = smooth(kgrid, disc, true); input_args = {&#39;PMLInside&#39;, false, &#39;PMLSize&#39;, PML_size, &#39;PlotPML&#39;, false, &#39;Smooth&#39;, false}; sensor_data_raw = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:}); % assign the time reversal data source2 = source; source2.p0=0; sensor2 = sensor; sensor2.time_reversal_boundary_data = sensor_data_raw; kgrid2 = kgrid; medium2 = medium; medium2.alpha_coeff = 0; p0_tr = kspaceFirstOrder2D(kgrid2, medium2, source2, sensor2, input_args{:});</code></pre>