k-Wave

1. SHickling
   

   Member
   Joined: May '13
   Posts: 9

   Hello,
   I am doing photoacoustic simulations with pulses that are ~5us, so the thermal confinement condition doesn't hold. I am assuming that the pulse is a rect function with constant amplitude throughout the entire pulse duration. I am trying to code this using source.p so each row corresponds to a source point and the columns are time steps (the amplitude of the pulse at each source location is different). The value that I am putting in for each entry is the total source pressure at that location divided by the number of time steps. My question though is how do I tell k-wave what the spacing for these time steps is?

   Thanks!
   Susannah

   Posted 11 years ago #

2. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,643

   Hi Susannah,

   The number and spacing for the time steps is set by kgrid.t_array. By default, after you call makeGrid this is set to 'auto'. In this case, the simulation functions then call makeTime to set the time array to something sensible. However, you can also set this yourself, for example using something like:

   kgrid.t_array = 0:dt:(Nt-1)*dt;

   where Nt is the desired number of time steps and dt is the spacing. There's some more information at the end of section 3.2 in the k-Wave Manual.

   Rather than using a rectangular input function, I'd recommend smoothing your input using filterTimeSeries, or using a smooth function like a Gaussian. This is because for a CFL < 1, your time discretisation will support more frequencies than your spatial discretisation, so using a temporal delta function (or anything with step-changes) will result in aliasing. There's some more information in the Numerical Analysis -> Filtering a Delta Function Input Signal Example in the MATLAB documentation.

   Brad.

   Posted 11 years ago #

3. SHickling
   

   Member
   Joined: May '13
   Posts: 9

   Hi Brad,

   Thanks for the response. I noticed on p. 36 of the k-Wave Manual that is says the source can have any number of time points and it doesn't need to be the same length as kgrid.t_array. I'm just confused about how I indicate what time step I am using to specify the pressure source, since it is different from t_array.

   Thanks,
   Susannah

   Posted 11 years ago #

4. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,643

   Hi Susannah,

   The time step for the pressure source and t_array must be the same, however, the number of time steps can vary. If the number of time points in the source is smaller than t_array, this means the source will stop transmitting before the end of the simulation.

   Brad.

   Posted 11 years ago #

5. SHickling
   

   Member
   Joined: May '13
   Posts: 9

   Hello again,

   I am seeing some strange behaviour when I try to use source.p and am wondering if you have any insight. I have an energy deposition pulse that lasts for 5us and have determined the pressure distribution that would result if the pulse was a delta function by using the Gruneisen coefficient. To account for the actual length of the pulse, I have tried to use source.p with the number of columns equal to (total pulse length)/(time_step), and each element being equal to the total amount of pressure divided by the number of time steps. However, the amplitude of the resulting acoustic pulse is about an order of magnitude larger when I do this compared to when the pressure is produced instantaneously. Does the source.p function do some sort of integration that I need to account for?

   Thanks!
   Susannah

   Posted 11 years ago #

6. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,643

   Hi Susannah,

   This is an interesting question. The difference is due to the scaling factor that is applied to time varying sources in k-Wave. This is outlined in Eqs. (2.19) and (2.17d) in the k-Wave manual. The total scaling factor applied is given by

   sf = [1 / (N * c0^2)] * [2 * c0 * dt / dx]

   where N is the number of dimensions, c0 is the sound speed, and dx and dt are the grid spacing and time step. For an initial pressure distribution (i.e., when using source.p0), only the first part of this scaling applies. This part of the scaling arises because within the code the pressure source terms are added to the acoustic density, where in the lossless case the relationship between pressure and density is given by:

   p = c0^2 * rho

   The N appears because the acoustic density is artificially divided into Cartesian components to allow a split-field perfectly matched layer (PML) to be implemented

   The second part of the scaling was added so that when using a time-varying plane wave source, the resulting amplitude in the simulation matches the input amplitude. To counter-act the second component of the scaling, you should divide your source magnitude by

   (2 * c0 * dt / dx)

   Also keep in mind, the longer you make your source, the less it will look like your instantaneous simulation.

   Hope that helps,

   Brad.

   Posted 11 years ago #

7. SHickling
   

   Member
   Joined: May '13
   Posts: 9

   That does the trick, thanks Brad!

   Posted 11 years ago #

8. kjones
   

   Member
   Joined: Oct '13
   Posts: 2

   Hello,

   I'm also trying to simulate the photoacoustic signal generated by longer pulses.
   I have run two simulations specifying source.p with different dt (but keeping sum(sum(source.p))*dt*dx*dy*dz constant). The amplitude of the results depends on dt, which it should not. Should source.p0 always be divided by (2 * c0 * dt / dx) to get "true" amplitudes?

   Kevin

   Posted 10 years ago #

9. bencox
   

   Developer
   Joined: Mar '10
   Posts: 292

   Hi Kevin,

   Excellent question. There is a reason for the behaviour you see.

   When your source is a single point in time, eg. source.p = 1, then the time series you measure at your detector will be the impulse response function. What this actually is will depend on the numerical method you are using. In k-Wave, for a homogeneous medium, you will see the convolution of the Green's function for the wave equation and the band-limited interpolant, see Section 2.8 in the manual. (What this looks like will depend on how many spatial dimensions you are working in.)

   When your source has two points in time, eg. source.p = [1 1], then what you will see will be the sum of two copies of the impulse response function spaced by dt in time. Usually dt is short enough that the two impulse responses will not be clearly separable, but will form an output with a single maximum, but at a higher level. (A useful exercise is to plot the measured time series as you increase the spacing between two pulses: source.p = [1 0 0 ... 0 0 1]. Eventually the two pulses measured at the detector will become separated in time.)

   Adding increasing numbers of time points to the source, source.p = [1 1 1 ... ], will continue to increase the maximum value of the output (well, there are some small oscillations) until eventually it settles down to a steady value. (Another useful exercise is to plot the maximum of the measured time series as a function of the number of pulses used in source.p.)

   So there are essentially two different regimes in which the 'source amplitude' is well-defined: a single pulse and a sufficiently long train of pulses. When using source.p, k-Wave by default scales the output so a plane wave excited by a long train of pulses will have the same amplitude as the source. When setting an initial value problem using source.p0 k-Wave does not include the scaling factor, so that, again, the output is in the same units as the initial pressure distribution.

   Spreading the source spatially can help ameliorate this effect (see Manual, Figure 2.6).

   Hope that helps,
   Ben

   Posted 10 years ago #

10. zahrat
    

    Member
    Joined: Mar '20
    Posts: 3

    Hello,

    I am having scaling issue as well and I am not sure if I understand your explanation correctly.
    I am simulating an ultrasound phased array, propagating along the z axis.

    My source is generated at follow:

    ampc = (5)*ones(size(phase));
    source.p = createCWSignals(t, freq, ampc,focDistRad, 0);%

    When I keep dx,dy=0.6 mm, sos=1485 fixed and just change dz for instance from .5mm to 0.25 mm, I am noticing the the output pressure max amp is reported double which I doubt it is due to the finer resolution.

    I get the dt from "maketime" fucntion (kgrid.makeTime(medium.sound_speed, 0.3)).

    You are mentioning that the scaling factor for the input (ie ampc) should be the inverse of (2 * c0 * dt / dz) while would nt it be the same number (ie .3) in my case (if I am using "maketime" function since dt= 0.3*sos/dz).

    Thanks for your help inadvance.
    Zahra

    Posted 4 years ago #

11. bencox
    

    Developer
    Joined: Mar '10
    Posts: 292

    Hi Zahra,

    The scaling factor of (2 * c0 * dt / dx) is exact for one-dimensional propagation, ie. a plane wave propagating in the x direction, but not for higher dimensions. In those cases it's a bit less straightforward. As we tried (!) to explain above, because of the numerical method that k-Wave is based on, when it takes an input 'source.p', it effectively multiplies each sample of that source with the band-limited interpolant, BLI - roughly sinc(pi*x/dx)sinc(pi*y/dy)sinc(pi*z/dz) - centred at the source point. Unlike the source point, the BLI is distributed in space, so when it propagates - which is what happens when a timestep is taken - there is a non-zero value of the field at the source point. The next source sample then adds to that existing field, and so on, and the result is that the effective amplitude of the wave generated by source.p is different from the input source.p. Furthermore, because the BLI depends on the grid spacing the amplitude depends on the grid spacing.

    One way to avoid this is to set source.p_mode = 'dirichlet', which enforces the values of source.p at the source points, but in this case the source points also act as scatterers while they are transmitting a signal, which can be a nuisance in some circumstances.

    A pragmatic alternative that I have used is to note how the amplitude scales with the grid size and incorporate that factor, eg. dz, into the expression for the source amplitude, so that as you change grid resolution the source amplitude remains constant.

    Hope that helps a bit.

    Best wishes
    Ben

    Posted 4 years ago #

12. zahrat
    

    Member
    Joined: Mar '20
    Posts: 3

    Hey Ben

    Thanks for your reply and clarification.

    I have two more questions. I am trying to compare the pressure value from the simulation with the measurement and they are not quite agreeing. I am reading lower with the simulation than expected. I am not confident that I set the input pressure (p0) correctly in the code. I have tried to input pressure and velocity and compare the results as well.

    1) One method I set the p0 as follow (power is in watt)

    %-----------------input pressure est
    pwr=2;Z0=1.5*1e6;area=1.35e-3*1.25e-3*6144;
    p0=sqrt(pwr*Z0/area); %---estimate of p0
    u0=p0/Z0; %----estimate of u0

    source.p_mask = tx_mask;
    ampc = (p0)*ones(size(phase));
    source.p_mode = 'dirichlet';

    source.p = createCWSignals(t, freq, ampc,focDistRad, 0);

    2) Using velocity approach as the source:

    source.u_mask = tx_mask;
    ampc_u=u0*ones(size(phase));
    source.u_mode = 'dirichlet';

    source.uz = createCWSignals(t, freq, ampc_u,focDistRad , 0);

    Two issues that I have is that 1) Max pressure from approach 1 & 2 are different 2) Max pressure from simulation seems to be lower than measurement.

    I appreciate if you let me know your thoughts on this.

    Thanks
    Zahra

    Posted 4 years ago #

13. bencox
    

    Developer
    Joined: Mar '10
    Posts: 292

    Hi Zahra,

    In 1D these sources should give the same result, but in higher dimensions they will differ as a pressure source is a monopole, which radiates omnidirectionally, but a velocity source is a dipole, which radiates with a figure-of-eight pattern. Could that be the reason for the difference?

    As for comparison to measurements, I would recommend having a look at our colleague Elly Martin's recent paper on comparing experiments with k-Wave. doi:10.1109/TUFFC.2019.2941795

    Best wishes
    Ben

    Posted 4 years ago #

14. HD Y
    

    Member
    Joined: Apr '23
    Posts: 6

    Hi, it seems I found some discussions on the implementation of Dirichlet boundary conditions in this post. But I still did not get how the scaling factor is set in general.

    On the other hand, I saw some good results in Figure 2 of the following publication:

    Modelling Nonlinear Ultrasound Propagation in Absorbing Media using the k-Wave Toolbox: Experimental Validation.

    May I ask in this reference, what exact scaling factor is used in order to get such good agreement between simulation results and experimental results.

    Thanks

    Posted 3 months ago #

RSS feed for this topic

Reply

You must log in to post.