k-Wave User Forum » Topic: compute green's function

k-Wave User Forum » Topic: compute green's function http://www.k-wave.org/forum/topic/compute-greens-function Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 00:17:31 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php bencox on "compute green's function" http://www.k-wave.org/forum/topic/compute-greens-function#post-6673 Wed, 05 Dec 2018 01:10:22 +0000 bencox 6673@http://www.k-wave.org/forum/ Hi Mohsen, Try turning the 'smooth' option off, by adding the optional argument <code>...,'Smooth', false)</code> in kspaceFirstOrder3D. Best wishes, Ben bencox on "compute green's function" http://www.k-wave.org/forum/topic/compute-greens-function#post-6672 Tue, 04 Dec 2018 13:34:56 +0000 bencox 6672@http://www.k-wave.org/forum/ Hi maryam.parto, Good question. I put that 4 in there to make the output equal 1, but now I look as it again it's not clear to me what that factor should be in the general case. I will investigate further.... Thanks for the question! Best wishes, Ben mohsen on "compute green's function" http://www.k-wave.org/forum/topic/compute-greens-function#post-6661 Mon, 26 Nov 2018 01:51:04 +0000 mohsen 6661@http://www.k-wave.org/forum/ Hello People to find a matrix model to build a linear system is in my interest. as you said I put the source.p0 = 1 which means an impulse in my pixel. for each pixel I put it and the solution which is impulse response will be assigned in the columns of the interested matrix. but when I multiply this matrix in the initial value pressure field although the behaviour of the plot is the same as the solution of k wave the magnitude is different. I will be thankful if somebody guides me the Best Mohsen maryam.parto on "compute green's function" http://www.k-wave.org/forum/topic/compute-greens-function#post-6573 Tue, 28 Aug 2018 10:17:16 +0000 maryam.parto 6573@http://www.k-wave.org/forum/ Hello dear bencox, My question is how do you scale the delta function for setting the initial pressure gradient? Actually, I can't understand the line "source.dp0dt=delta./4dt". why did you divide it by 4? Generally, I would be grateful if you explain me how can I obtain the initial pressure gradient from a known initial pressure. thanks a lot bencox on "compute green's function" http://www.k-wave.org/forum/topic/compute-greens-function#post-342 Thu, 01 Mar 2012 00:56:37 +0000 bencox 342@http://www.k-wave.org/forum/ Hi Chao, Yes, very similar. Both will remove the high frequency components. Ben huangchao on "compute green's function" http://www.k-wave.org/forum/topic/compute-greens-function#post-341 Thu, 01 Mar 2012 00:45:08 +0000 huangchao 341@http://www.k-wave.org/forum/ Hi, Dr. Cox, I have read your paper, and it's clearer to me now about the oscillation. In your paper, the oscillations in the recorded time series were reduced by smoothing the initial pressure distribution (Kronecker delta function) by using windows in spatial frequency domain (k-space), and I am wondering if this is equivalent to smooth the recorded time series by use of windows in temporal frequency domain. Thanks, Chao bencox on "compute green's function" http://www.k-wave.org/forum/topic/compute-greens-function#post-330 Fri, 24 Feb 2012 10:21:37 +0000 bencox 330@http://www.k-wave.org/forum/ Hi Chao, The oscillations come from the fact that what you are seeing is the propagation of the underlying bandlimited function, which (as the code in my previous post showed) is an oscillating function. So you are getting what you should expect to get. It is not that it is inaccurate, it is that what you thought was originally a delta function wasn't (as explained above) it was just that at t=0 you were sampling it mostly at zeros so it looked like a delta function. At other times you might be sampling at different points on the wave form and so you will now see oscillations. Have a look at <a href="http://www.k-wave.org/papers/2011-Treeby-JASA.pdf">this paper</a> for more explanation. As I said before, you can only get a bandlimited version of the Green's function when calculating it discretely. The only way to get closer to a Dirac delta function is to user a finer mesh. If you define your discrete delta properly (so its height grows as the pixel size reduces) then in the limit of an infinitesimally small pixel size it will converge towards the Dirac delta function, and the output will converge towards the Green's function (without any bandlimiting). Hope that's clearer, Ben huangchao on "compute green's function" http://www.k-wave.org/forum/topic/compute-greens-function#post-329 Thu, 23 Feb 2012 22:06:53 +0000 huangchao 329@http://www.k-wave.org/forum/ Hi Dr. Cox, Thank you very much for the detailed explanation. It's very helpful. I have one more question, that is when I set source.p0 to be 1 at a single grid point to calculate the impulse response ('sensor_data' in the code), or the time derivative of the Green's function, whatever you call it, it has a lot of oscillations. Of course, this can be avoid by smoothing the source.p0, but it will compromise the accuracy of the computed impulse response, because the source is not Kronecker δ after smoothing. Any ideas to how to avoid the oscillation and also maintain the impulse response accuracy at the same time? Thanks, Chao bencox on "compute green's function" http://www.k-wave.org/forum/topic/compute-greens-function#post-328 Thu, 23 Feb 2012 12:45:54 +0000 bencox 328@http://www.k-wave.org/forum/ Hi Chao, On point (1): what you think is a perfect delta function, because you only see a single 1 and many zeros, can be thought of as a discrete representation of a continuous signal, i.e. you are only looking at a sampled version of the underlying continuous signal. The underlying continuous signal depends on the basis functions you choose (i.e. how you choose to interpolate between the samples) and for a Fourier basis, the underlying function that fits these samples is a bandlimited, oscillatory function. Run this code and you will see what I mean. <code> % create the computational grid Nx = 64; % number of grid points in the x (row) direction dx = 1; % grid point spacing in the x direction [m] x = (-Nx/2:Nx/2-1)*dx; % define a single source pixel (numerical delta function) delta = zeros(Nx,1); delta(Nx/2+1) = 1; % calculate underlying bandlimited delta function using a much finer mesh dx_fine = 0.001; x_fine = min(x):dx_fine:max(x); underlying_func = (1/Nx)*cot(x_fine*pi/(dx*Nx)).*sin(x_fine*pi/dx); underlying_func(x_fine==0) = 1; % plots figure subplot(3,1,1) stem(x,delta,'o') axis([min(x) max(x) -0.22 1]) title('''discrete'' delta function') subplot(3,1,2) plot(x_fine,underlying_func,'k') axis([min(x) max(x) -0.22 1]) title('underlying bandlimited function') subplot(3,1,3) stem(x,delta,'o') hold on plot(x_fine,underlying_func,'k') axis([min(x) max(x) -0.22 1]) title('''discrete'' delta is sampled version of bandlimited function') </code> (2) Yes. You could think of this as a linear system (as long as you stick to the linear wave equation and not the nonlinear version). You can imagine a matrix equation with a vector of p(x,t) stacked up on the left, a matrix made up by discretising the wave operator, and an initial condition (or a source, depending on how you have written the operator matrix) on the right. You haven't avoided solving the wave equation: this just is a method for solving the wave equation (and probably not a very efficient one). You need to choose some way of discretizing the wave operator to form the system matrix, and that might be using finite differences, finite elements, spectral approaches, or whatever you like! (3) These are equivalent in the sense that p_0 = \Gamma H is true for photoacoustics (when the optical pulse is short enough), and dp/dt = 0 when u = 0 (which can be seen by combining the mass conservation equation and pressure-density relation). Hope that helps, Ben huangchao on "compute green's function" http://www.k-wave.org/forum/topic/compute-greens-function#post-327 Thu, 23 Feb 2012 07:02:39 +0000 huangchao 327@http://www.k-wave.org/forum/ Dear Dr. Cox, Thanks for your reply, and I have some additional questions: 1. I don't understand why Kronecker delta is a band-limited version of the Dirac delta function. In my opinion, the spectrum of Kronecker delta is still infinitely wide. For example, if you plot the amplitude spectrum of the Kronecker delta using MATALB commands 'plot(abs(fft([1, zeros(1,N)])))', the amplitude of all frequency components is 1. 2. I agree that setting source.p0 is 1 at a single grid point will give the time derivative of the Green's function. This leads to another relevant question: Is solving the linear wave equation (no matter the 2nd order equation, or the 3 coupled 1st order equations) is equivalent to compute the output of a linear system where source is the input? If so, for the photoacoustic wave equation: <img src="http://i.imgur.com/C3jyF.jpg" /> we can treat the source H as the input, p as the output, and other parts of the equation (including the time derivative on the right side) as the linear system operator. In this case, because H=A(r)δ(t) is used as the source term for photoacoustic wave equation, just setting source.p0 is 1 at a single grid point (A(r) becomes Kronecker δ) is equivalent to setting the input as an impulse (if Kronecker δ is indeed a good approximation of Dirac δ). Then using k-Wave to solve the equation will give us the impulse response of the linear system, which can be used to compute the output (p) just like any other linear systems if somehow we don't want to solve the wave equation to compute p. 3. I noticed that when solving the wave equation, 2nd order equation or 1st order equations, there are 2 kinds of initial conditions: <img src="http://i.imgur.com/ZVpQ9.jpg" /> <img src="http://i.imgur.com/fvZT5.jpg" /> Could you tell me what's the relation between these 2 initial conditions? Are they equivalent? Thanks again for your help! Best, Chao bencox on "compute green's function" http://www.k-wave.org/forum/topic/compute-greens-function#post-325 Wed, 22 Feb 2012 14:02:37 +0000 bencox 325@http://www.k-wave.org/forum/ Dear Chao, The frequency spectrum of a Dirac delta function is infinitely wide, whereas any numerical model will have some upper limit on the frequencies it can deal with which is set by the spacing of the grid, so the best you can hope for is a band-limited version of the Green's function. (When you input a Kronecker delta you are effectively using a band-limited version of the Dirac delta function). In practice, you just need to ensure the bandwidth is wide enough for your application. Setting the initial condition that <code>source.p0</code> is 1 at a single pixel will not give you the Green's function but the time derivative of the Green's function, as that initial condition is equivalent to using the time derivative of the delta function as a source term. The Green's function is defined as the solution that arises from using a delta function as the source term which you can get by setting the initial pressure gradient. Something like this example below (which easily extends to 2D and 3D) gives the gist of it. <code> % create the computational grid Nx = 512; % number of grid points in the x (row) direction dx = 0.1e-3; % grid point spacing in the x direction [m] kgrid = makeGrid(Nx, dx); k = kgrid.k; % define the properties of the propagation medium medium.sound_speed = 1500; % [m/s] % define a time array [kgrid.t_array dt] = makeTime(kgrid,medium.sound_speed); kgrid.t_array = (0:length(kgrid.t_array)/2)*dt; % define a single source pixel delta = zeros(Nx,1); delta(Nx/2+1) = 1; %scale the delta function to have unit area (in time) and assign to the %initial pressure gradient source.dp0dt = delta./(4*dt); % define a single sensor point sensor.mask = zeros(Nx,1); sensor.mask(Nx/4) = 1; % run the simulation output = kspaceSecondOrder(kgrid, medium, source, sensor,'PlotScale',[0 1.5]); %plot the output, which is the bandlimited Green's function. figure plot(kgrid.t_array, output) </code> On your other point, k-Wave does the spatial calculations using Fourier basis functions (it is a type of Fourier pseudospectral method), but there won't be any additional higher frequency information between the pixels in the grid so you could just interpolate linearly and it will be ok. Ben huangchao on "compute green's function" http://www.k-wave.org/forum/topic/compute-greens-function#post-323 Wed, 22 Feb 2012 07:56:48 +0000 huangchao 323@http://www.k-wave.org/forum/ Hi, I want to compute the wave filed from a point source, which is called Green's function or point spread function, by using k-Wave. Theoretically, the point source is described by a Dirac delta function, so I don't know if it's appropriate to use Kronecker delta function to approximate Green's function in k-Wave to compute Green's function. Using Kronecker delta function, I mean to set 1 pixel value of the source.p0 to be 1, and other pixel values to be zeros. Another relevant question is what the basis or expansion function is used in k-Wave. I mean if I have the initial pressure distribution source.p0 as a matrix, but I want to know the pressure between grid points, how to compute that from source.p0? Thanks in advanced! Best, Chao