http://www.k-wave.org/forum/topic/diagonal-medium-frequency-artifact-in-helmholtz-inversion Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 14:17:04 +0000 http://bbpress.org/?v=1.0.2 q http://www.k-wave.org/forum/search.php http://www.k-wave.org/forum/topic/diagonal-medium-frequency-artifact-in-helmholtz-inversion#post-5937 Sun, 21 May 2017 16:16:45 +0000 Bradley Treeby 5937@http://www.k-wave.org/forum/ <p>Hi Eric,</p> <p>I have run your code, but to be honest, I'm not really sure what I'm looking for. Do you suspect it is something wrong with the wave simulation, or the post-processing? If it's the latter, I probably won't be much help I'm afraid, as I don't know a lot about the specific problem you're working on.</p> <p>Brad. </p> http://www.k-wave.org/forum/topic/diagonal-medium-frequency-artifact-in-helmholtz-inversion#post-5837 Tue, 21 Feb 2017 17:34:40 +0000 Dr. Eric Barnhill 5837@http://www.k-wave.org/forum/ <p>Hi Bradley, thank you for your interest in this problem. I have tested the following matlab script and it should run needing only kwave on the path. It inlines all the necessary functions, and replicates the condition I stated. You can run it by calling</p> <p>G = kwave_inversion_demo;</p> <p>then inspecting the object G which is the shear modulus. The correct answer is 4000 and it is mostly very close, however you will see some smooth artifact comes through that is far off, and needless to say, the material is supposed to be homogeneous.</p> <p>Okay here is the code:<br /> -----<br /> <pre><code>function G = kwave_inversion_demo K = 5; G = 2; dx = .002; f = 50; ncycles = 4&#39;; prefs = kwave_prefs(&#39;K&#39;, K, &#39;G&#39;, G, &#39;dx&#39;, dx, &#39;f&#39;, f, &#39;ncycles&#39;, ncycles); u = kwave_eb(prefs); u = cat(5, u.ux, u.uy, u.uz); u_t = size(u,4); final_cycle_start = round(u_t*(ncycles-1)/ncycles); u_cycle = u(:,:,:,final_cycle_start:end,:); % temporal fourier transform, retain driving frequency bin u_ft = fft(u_cycle, [], 4); u_ft = squeeze(u_ft(:,:,:,2,:)); % apply Helmholtz equation to complex wavefield u_lap = lap(u_ft) ./ (prefs.dx)^2; G = prefs.rho*(2*pi*prefs.f)^2*sum(abs(u_ft), 4) ./ sum(abs(u_lap), 4); end function u = kwave_eb(prefs) SIZE = 96; PML_SIZE = 16; SOURCE_POS = SIZE / 2; Nx = SIZE; % number of grid points in the x direction Ny = SIZE; % number of grid points in the y direction Nz = SIZE; % number of grid points in the z direction dx = prefs.dx; % grid point spacing in the x direction [m] dy = prefs.dx; % grid point spacing in the y direction [m] dz = prefs.dx; % grid point spacing in the z direction [m] kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz); medium.sound_speed_compression = prefs.K*ones(Nx, Ny, Nz); % [m/s] medium.sound_speed_shear = prefs.G*ones(Nx, Ny, Nz); % [m/s] medium.density = prefs.rho*ones(Nx, Ny, Nz); % [kg/m^3] % define the absorption properties medium.alpha_coeff_compression = 0.1; % [dB/(MHz^2 cm)] medium.alpha_coeff_shear = 0.1; % [dB/(MHz^2 cm)] % create the time array t_end = prefs.t; kgrid.t_array = makeTime(kgrid, max(medium.sound_speed_compression(:)), [], t_end); source.u_mask = zeros(Nx, Ny, Nz); source.u_mask(SOURCE_POS, SOURCE_POS, SOURCE_POS) = 1; % define source to be a velocity source source_freq = prefs.f; % [Hz] source_mag = 1; % [Pa] source.ux = source_mag*sin(2*pi*source_freq*kgrid.t_array); source.uy = source_mag*sin(2*pi*source_freq*kgrid.t_array); source.uz = source_mag*sin(2*pi*source_freq*kgrid.t_array); %source.ux = filterTimeSeries(kgrid, medium, source.ux); %source.uy = filterTimeSeries(kgrid, medium, source.uy); %source.uz = filterTimeSeries(kgrid, medium, source.uz); % define sensor mask in x-y plane using cuboid corners, where a rectangular % mask is defined using the xyz coordinates of two opposing corners in the % form [x1, y1, z1, x2, y2, z2].&#39; sensor.mask = [1 + PML_SIZE, 1 + PML_SIZE, 1 + PML_SIZE, Nx - PML_SIZE, Ny - PML_SIZE, Nz - PML_SIZE].&#39;; % record the maximum pressure in the plane sensor.record = {&#39;u&#39;}; % define input arguments input_args = {&#39;PMLSize&#39;, PML_SIZE}; % run the simulation with PML inside u = pstdElastic3D(kgrid, medium, source, sensor, input_args{:}); end function prefs = kwave_prefs(varargin) prefs = default_prefs; fields = get_fields; for n = 1:2:nargin name = varargin{n}; value = varargin{n+1}; for f = 1:numel(fields) if strcmpi(fields{f}, name) evalc([&#39;prefs.&#39;,name,&#39;=&#39;,num2str(value),&#39;;&#39;]); end end end end function prefs = default_prefs prefs.K = 1540; prefs.G = 1.54; prefs.dx = 0.002; prefs.rho = 1000; prefs.f = 50; prefs.ncycles = 4; prefs.t = (1/prefs.f)*prefs.ncycles; end function fields = get_fields fields = {&#39;K&#39;, &#39;G&#39;, &#39;dx&#39;, &#39;rho&#39;, &#39;f&#39;, &#39;ncycles&#39;, &#39;t&#39;}; end function y = lap(x) laplacian = zeros(3,3,3); laplacian(2,2,1) = 1; laplacian(2,2,3) = 1; laplacian(:,:,2) = [0 1 0; 1 -6 1; 0 1 0]; y = convn(x, laplacian, &#39;same&#39;); end</code></pre> http://www.k-wave.org/forum/topic/diagonal-medium-frequency-artifact-in-helmholtz-inversion#post-5832 Thu, 16 Feb 2017 23:09:11 +0000 Bradley Treeby 5832@http://www.k-wave.org/forum/ <p>Hi Eric,</p> <p>It's hard to say without seeing your code. Could you post an example that recreates the artefact that you mention?</p> <p>Brad. </p> http://www.k-wave.org/forum/topic/diagonal-medium-frequency-artifact-in-helmholtz-inversion#post-5830 Wed, 15 Feb 2017 15:53:58 +0000 Dr. Eric Barnhill 5830@http://www.k-wave.org/forum/ <p>First of all thanks for the excellent software. My first intended use for this toolbox is shear wave inversion assuming homogeneous linear elasticity. So, in its simplest form the field can just be divided by its Laplacian, multiplied by some constants. </p> <p>However I get a funny artifact when I do this, a medium-spatial frequency artifact that goes from one diagonal of the simulation to the other. It is like an oscillating pie wedge. I have posted a picture of it at <a href="http://imgur.com/a/K4NrA" rel="nofollow">http://imgur.com/a/K4NrA</a>. </p> <p>The compression modulus is set &gt;20x the shear, and neglected. Could that be the source of this artifact? </p>