http://www.k-wave.org/forum/topic/simple-3d-simulation-received-signal-question Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 00:47:51 +0000 http://bbpress.org/?v=1.0.2 q http://www.k-wave.org/forum/search.php http://www.k-wave.org/forum/topic/simple-3d-simulation-received-signal-question#post-4759 Wed, 08 Oct 2014 21:48:15 +0000 Bradley Treeby 4759@http://www.k-wave.org/forum/ <p>Hi Margot,</p> <p>I see exactly what you mean, but I'm afraid I'm no closer to an explanation! Would you be happy to share (or email me) your code?</p> <p>Brad. </p> http://www.k-wave.org/forum/topic/simple-3d-simulation-received-signal-question#post-4758 Mon, 06 Oct 2014 09:07:03 +0000 mgt_p 4758@http://www.k-wave.org/forum/ <p>Hi Brad,</p> <p>My transducer.focus_distance is .04 + dx (so .04025m). I have not set any transducer.elevation_focus_distance (it makes it Inf, right?), and no steering. </p> <p>I plotted the raw sensor_data with the same scheme as before (I cut the signal in two, and normalize each half independantly so that the input signal and the scatterer echo both appear with an amplitude of 1). Here it is:<br /> <a href="https://www.dropbox.com/s/cydrfgx5b87wwye/not_preprocessed_sensor_data.PNG?dl=0" rel="nofollow">https://www.dropbox.com/s/cydrfgx5b87wwye/not_preprocessed_sensor_data.PNG?dl=0</a><br /> <a href="https://www.dropbox.com/s/u0a71s4tn9rmks6/raw_sensor_data.fig?dl=0" rel="nofollow">https://www.dropbox.com/s/u0a71s4tn9rmks6/raw_sensor_data.fig?dl=0</a></p> <p>It shows quite clearly that the input signal has the drift already... </p> <p>Regards,</p> <p>Margot </p> http://www.k-wave.org/forum/topic/simple-3d-simulation-received-signal-question#post-4754 Fri, 03 Oct 2014 16:35:30 +0000 Bradley Treeby 4754@http://www.k-wave.org/forum/ <p>Hi Margot,</p> <p>Could you tell us a bit more about the exact simulation setup? What are your settings for <code>transducer.focus_distance</code>, <code>transducer.elevation_focus_distance</code>, and <code>transducer.steering_angle</code>? Does the raw <code>sensor_data</code> look ok before processing? </p> <p>You can include images in your posts as explained <a href="http://www.k-wave.org/forum/topic/using-images-in-posts">here</a>, but you need to host them somewhere else I'm afraid.</p> <p>Brad. </p> http://www.k-wave.org/forum/topic/simple-3d-simulation-received-signal-question#post-4746 Fri, 26 Sep 2014 16:46:36 +0000 mgt_p 4746@http://www.k-wave.org/forum/ <p>Hi Brad,</p> <p>New simulation, new problem.</p> <p>I am now using an array that is more realistic than the previous one, with 64 elements and plausible size, but the rest is more or less the same.</p> <p>I am looking at sensor_data, just after the simulation, since I want to investigate receive beamforming.</p> <p>I can observe a drift (although not very linear) along my channels, as in my first echo appears later in sensor_data(64,:) than in sensor_data(1,:) when they have exactly the same spatial/physical phenomenon origin.</p> <p>You can have a look at plots of side-by-side channels here:<br /> <a href="https://www.dropbox.com/s/827cd4yk7fikd9r/side_by_side_sensor_data.PNG?dl=0" rel="nofollow">https://www.dropbox.com/s/827cd4yk7fikd9r/side_by_side_sensor_data.PNG?dl=0</a><br /> (all amplitudes should be disregarded, I do a lot of normalizations here and there to ease the visualization)<br /> And if I zoom on the interface (which is the layer with higher sound speed), here is what I see:<br /> <a href="https://www.dropbox.com/s/kyh610hpqh7j64x/interface_drift.PNG?dl=0" rel="nofollow">https://www.dropbox.com/s/kyh610hpqh7j64x/interface_drift.PNG?dl=0</a></p> <p>However, appart from the normalization, I do get rid of the input signal using the Tukey windowing procedure used in the k-wave examples for medical ultrasound, and crop them to spatially relevant zones (from 0m to the actual end of my phantom).</p> <p>I can't figure out an explanation for this drift, my hypotheses being:<br /> - the sensors do not start "recording" at the same time<br /> - something happens in transmit beamforming (but I haven't looked into how it's done in k-wave at all).</p> <p>Do you have any hint as to what might be happening here? Or a direction in which I should look? I absolutely need to get rid of that drift (and without interpolation if possible)...</p> <p>Thank you very much,</p> <p>Margot</p> <p>PS: Is there a way to include the images in the posts in the forum? </p> http://www.k-wave.org/forum/topic/simple-3d-simulation-received-signal-question#post-4686 Wed, 20 Aug 2014 07:52:14 +0000 mgt_p 4686@http://www.k-wave.org/forum/ <p>Hi Brad,</p> <p>Thanks for the additional information.<br /> Yes, I had tried the C++ code when I was running bigger simulations, but it seemed kind of less reliable to work with non-beamformed data (cf my other post in the forum), so this time I decided I could afford the time :D But I will, eventually, go back to it I think.</p> <p>Regards,</p> <p>Margot </p> http://www.k-wave.org/forum/topic/simple-3d-simulation-received-signal-question#post-4684 Tue, 19 Aug 2014 17:59:24 +0000 Bradley Treeby 4684@http://www.k-wave.org/forum/ <p>Hi Margot,</p> <p>This factorisation is fine. However, the small prime factors should be for the <em>total</em> grid size, including the PML if <code>PMLInside</code> is set to <code>false</code>. Regarding the compute times, have you tried using the C++ code? It is relatively simple to use, you just need to download the binary and replace your calls to <code>kspaceFirstOrder3D</code> with <code>kspaceFirstOrder3DC</code>.</p> <p>Brad. </p> http://www.k-wave.org/forum/topic/simple-3d-simulation-received-signal-question#post-4683 Tue, 19 Aug 2014 16:05:34 +0000 mgt_p 4683@http://www.k-wave.org/forum/ <p>Dear Brad,</p> <p>Thank you for your answer, as a matter of fact, I was also on holidays, so the timing did not matter for me!</p> <p>About the values of Nx, Ny and Nz, I have some reality constraints that I can't play too much with. Nevertheless, I tried to have better prime factors (2 2 2 2 2 7 for example), but I still end up with over 1hr of simulation time. I wonder then if the PML should be added for the prime factor effect to appear? Or is this decomposition not small enough?</p> <p>Thank you,</p> <p>Margot </p> http://www.k-wave.org/forum/topic/simple-3d-simulation-received-signal-question#post-4680 Mon, 18 Aug 2014 18:17:46 +0000 Bradley Treeby 4680@http://www.k-wave.org/forum/ <p>Hi Margot,</p> <p>Apologies for taking so long to reply - it has been summer holidays here. </p> <p>The issue is related to high frequencies in your tone burst not being supported by the grid. You could resolve this by lowering the source frequency, decreasing the grid spacing, or filtering out the high frequencies from your input signal. </p> <p>For example, to do the latter you could add the following code after scaling:</p> <pre><code>% add some zeros to the start and end of the signal to give the % filter room to work input_signal = [zeros(size(input_signal)), input_signal, zeros(size(input_signal))]; % plot the amplitude spectrum before filtering spect(input_signal, 1/kgrid.dt, &#39;Plot&#39;, [true, false]); % remove high frequency components med.sound_speed = 1500; input_signal = filterTimeSeries(kgrid, med, input_signal, &#39;ZeroPhase&#39;, true); % plot the amplitude spectrum after filtering spect(input_signal, 1/kgrid.dt, &#39;Plot&#39;, [true, false]);</code></pre> <p>Adding this gave me nice clean looking reflections when plotting <code>sensor_data(5, :)</code> after running your code snippet. </p> <p>I'd advise using <code>filterTimeSeries</code> for all time varying input signals to be on the safe side. Notice I added some zeros before and after the signal to give it room to work (removing the high frequencies increases the width of the signal). I have also forced it to be zero phase so the features are not shifted as they would be if the filter was causal.</p> <p>For your simulation, I'd also suggest setting the values of Nx, Ny, and Nz explicitly so that they have small prime factors (you can check this using <code>checkFactors</code> or <code>factor</code>). This will make the simulations much faster.</p> <p>Hope that helps,</p> <p>Brad. </p> http://www.k-wave.org/forum/topic/simple-3d-simulation-received-signal-question#post-4657 Wed, 23 Jul 2014 12:59:27 +0000 mgt_p 4657@http://www.k-wave.org/forum/ <p>Dear Mr Treeby,</p> <p>I am trying to simulate a "simple" set up, and I don't really understand the signals I receive.</p> <p>My phantom is homogenous in density but for a 2x2x2 pixels highly reflective scatterer in the middle (high density). In terms of speed of sound, it is homogenous, with a bulk at 1540m/s and a layer, placed before the scatterer, at higher speed of sound. The scatterer gets an extra reflective punch from having also a higher speed of sound. There is no attenuation set.</p> <p>I have a linear transducer, with only 10 elements, focused on my scatterer in azimuth, and not focused in elevation.</p> <p>My first question is regarding the signal I get when I do not inforce the layer.<br /> I get this:<br /> <a href="https://www.dropbox.com/s/cwymc0crk3tn9cc/Signal_no_layer_after_Tukey_before_TGC.PNG" rel="nofollow">https://www.dropbox.com/s/cwymc0crk3tn9cc/Signal_no_layer_after_Tukey_before_TGC.PNG</a></p> <p>This, as the name says, is without any other post-processing step than the Tukey windowing to remove the trace of the input signal, but cropped to a the depth of my phantom.<br /> I do not understand the first, very big, echo, as my scatterer shows up where it is expected, at 2cm*. Do you know where this comes from?</p> <p>* about the scatterer position: I actually found out that if I wanted it to be at the expected depth, I had to build up the r vector not with t0 at half the input signal, but at a quarter of it. Any idea why? I had the feeling that the half-signal correction was more like a compromise than a definite rule.</p> <p>Then, when I put up the layer, I see echos within it that I'm not expecting.<br /> <a href="https://www.dropbox.com/s/mb8o57iayypaaty/Signal_layer_after_Tukey_before_TGC.PNG" rel="nofollow">https://www.dropbox.com/s/mb8o57iayypaaty/Signal_layer_after_Tukey_before_TGC.PNG</a><br /> The fact that my scatterer is drowned in the tail is not so much of a problem if I can explain where this "tail" comes from. But the presence of some signal between the two interfaces of the layer is strange to me. </p> <p>My hypothesis at the moment is that the grid spacing is too coarse. Do you have a better idea?</p> <p>Thank you very much,</p> <p>Margot</p> <p>Code:<br /> <pre><code>% define heterogenous medium % phantom size X = 0.04; %[m] Y = 0.02; %[m] Z = 0.02; %[m] % SOS bulk c = 1540; % [m/s] delta_c = 18; % [m/s] % density bulk rho = 1000; % [kg/m^3] % highly reflective scatterer properties scatterer.c = 1590; scatterer.rho = 2000; % [kg/m^3] % grid spacing dx = c/(3*2.85e6); % [m] dy = dx; % better for stability dz = dx; % better for stability % define the transducer % physical properties transducer.number_elements = 10; transducer.element_width = round(9.905e-4/dx); transducer.element_length = transducer.element_width; transducer.element_spacing = 0; transducer.radius = Inf; % effective grid size Nx = round(X/dx); Nz = round(Z/dz); n_scan_lines = 1; Ny_tot = round(Y/dy); Ny = Ny_tot - (n_scan_lines-1)*transducer.element_width; % layer depth L_top = 0.01; %[m] LTH = 0.005; %[m] layer thickness % highly reflective scatterer position scatterer.x = X/2; %[m] scatterer.y = Y/2; %[m] scatterer.z = Z/2; %[m] % PML size for border conditions PML_X = 30; PML_Y = 15; PML_Z = 15; % make grid % space kgrid = makeGrid(Nx,dx,Ny,dy,Nz,dz); % time t_end = (Nx*dx)*2.2/c; [kgrid.t_array,dt] = makeTime(kgrid, c, [], t_end); % medium physical properties SOS_map = repmat(c,[Nx,Ny_tot,Nz]); rho_map = repmat(rho,[Nx,Ny_tot,Nz]); % layer mask layer_mask = zeros(size(SOS_map)); layer_mask(round(L_top/dx):round(L_top/dx)+round(LTH/dx),:,:) = 1; % neutral background map SOS_map(layer_mask == 1) = SOS_map(layer_mask == 1)+ delta_c; % highly reflective scatterer SOS_map(round(scatterer.x/dx):round(scatterer.x/dx)+1,... round(scatterer.y/dy):round(scatterer.y/dy)+1,... round(scatterer.z/dz):round(scatterer.z/dz)+1) = scatterer.c; rho_map(round(scatterer.x/dx):round(scatterer.x/dx)+1,... round(scatterer.y/dy):round(scatterer.y/dy)+1,... round(scatterer.z/dz):round(scatterer.z/dz)+1) = scatterer.rho; % Transducer geometry % position and size element_size = [transducer.element_width transducer.element_length]*dx; transducer_aperture = element_size(1)*transducer.number_elements + (transducer.number_elements-1)*transducer.element_spacing*dx; transducer_width = transducer.number_elements*transducer.element_width ... + (transducer.number_elements - 1)*transducer.element_spacing; transducer.position = round([1, Ny/2 - transducer_width/2, Nz/2 - transducer.element_length/2]); % beamforming properties transducer.sound_speed = c; transducer.focus_distance = X/2 + dx; % in the middle of the phantom transducer.steering_angle = 0; transducer.active_elements = ones(transducer.number_elements,1); transducer.transmit_apodization = &#39;Hanning&#39;; transducer.receive_apodization = &#39;Hanning&#39;; % input signal % define properties of the input signal source_strength = 1e6; % [Pa] tone_burst_freq = 2.85e6; % [Hz] tone_burst_cycles = 4; % create the input signal using toneBurst input_signal = toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst_cycles); % scale the source magnitude by the source_strength divided by the % impedance (the source is assigned to the particle velocity) input_signal = (source_strength./(c*rho)).*input_signal; % assign transducer.input_signal = input_signal; % create transducer = makeTransducer(kgrid,transducer); % Simulation scan scan_line = zeros(n_scan_lines,kgrid.Nt); medium_position = 1; for l = 1:n_scan_lines medium.sound_speed = SOS_map(:,medium_position:medium_position + Ny - 1 ,:); medium.density = rho_map(:,medium_position:medium_position + Ny - 1,:); input_args = {... &#39;PMLInside&#39;, false, &#39;PlotSim&#39;, false, &#39;PMLSize&#39;, [PML_X, PML_Y, PML_Z], ... &#39;DataCast&#39;, &#39;single&#39;}; sensor_data = kspaceFirstOrder3D(kgrid,medium,transducer,transducer,input_args{:}); [rf,ch_data] = transducer.scan_line(sensor_data);</code></pre>