k-Wave

1. huangchao
   

   Member
   Joined: Nov '10
   Posts: 57

   Hi,

   I want simulate a plane wave, is there any other ways to do this besides using a plane source (multiple point sources on a very large plane)?

   Thank you very much!

   Best,
   Chao

   Posted 13 years ago #

2. bencox
   

   Developer
   Joined: Mar '10
   Posts: 292

   Hi Chao,
   At the moment that is the way to do it. If the whole problem is planar then you could do it in 1D, which would be much lighter on the memory.
   Ben

   Posted 13 years ago #

3. huangchao
   

   Member
   Joined: Nov '10
   Posts: 57

   Hi Dr. Cox,

   Thank you for your suggestion. I think simulating plane wave in 1D is a good idea, but I found that the signal will decay in 1D simulation (no absorption), as shown below

   

   I don't why the amplitude will decay in 1D if it's plane wave.

   Thanks,
   Chao

   Posted 13 years ago #

4. bencox
   

   Developer
   Joined: Mar '10
   Posts: 292

   Hi Chao,
   It shouldn't do, if the absorption's turned off. Can you post your code?
   Thanks,
   Ben

   Posted 13 years ago #

5. huangchao
   

   Member
   Joined: Nov '10
   Posts: 57

   Hi Dr. Cox,

   Just ignore my last post. The amplitude does not decay. It's my mistake. Sorry about that.

   Chao

   Posted 13 years ago #

6. huangchao
   

   Member
   Joined: Nov '10
   Posts: 57

   Hi Dr. Cox,

   Sorry to bother you again, but I noticed a phenomenon in 1D simulation with heterogeneous medium, and I don't know if it's reasonable.

   Nx = 512;
   dx = 0.1e-3;
   kgrid = makeGrid(Nx, dx);

   medium.sound_speed = 1500*ones(Nx,1);
   medium.sound_speed(200:300) = 3000;
   medium.density = 1000*ones(Nx,1);
   medium.density(200:300) = 2000;

   [kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed);
   kgrid.t_array = (0:4999)*dt;

   x_pos = 70;
   width = 20;
   height = 1;
   in = 0:pi/(width/2):2*pi;
   source.p0 = [zeros(x_pos,1); ((height/2)*sin(in-pi/2)+(height/2)).'; zeros(Nx - x_pos - width - 1, 1)];

   sensor.mask = zeros(Nx, 1);
   sensor.mask(150) = 1;
   sensor.mask(400) = 1;

   sensor_data = kspaceFirstOrder1D(kgrid, medium, source, sensor, 'PlotSim', false, 'Smooth', false);

   The signal from the first sensor looks like:
   
   The signal from the second sensor looks like:
   

   I don't know why the reflect wave pressure recorded in the first sensor become negative, while the signals from the 2nd sensor are always positive. Is this reasonable?

   Thank you very much!

   Best,
   Chao

   Posted 13 years ago #

7. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,643

   Hi Chao,

   When a wave propagates from one medium to another, if the second medium has a higher impedance, the reflected wave will be in phase with the incident wave. If the second medium has a lower impedance, the reflected wave will instead be 180 degrees out of phase (so a positive pressure wave will become a negative pressure wave and vice-versa). In both cases, the transmitted wave will be in phase with the incident wave.

   In your example, you have a region of higher impedance in the middle of your medium. In this case, the wave that is reflected back and forth will change signs each time it reaches the medium interfaces because the surrounding medium has a lower impedance. You can see more clearly what is happening if you switch the plot setting on.

   I hope that helps. You can find more information in a standard acoustics text, for example Kinsler & Frey (Page 152 in the 4th Edition).

   Brad.

   Posted 13 years ago #

8. rarslana
   

   Member
   Joined: Jun '12
   Posts: 4

   Hi,

   I am trying to define two sources with different frequencies propagating in the same kgrid. In order to do this I have defined two source.p_mask as follow
   source.p_mask = zeros(Nx, Ny);
   source.p_mask(end - Nx*0.625, Ny*0.375) = 1;
   source.p_mask(end - Nx*0.375, Ny*0.375) = 1;

   and in order to assign sources with different frequency to each of this source mask I did as follow;

   source.p(end - Nx*0.625, Ny*0.375) = source_mag*sin(2*pi*source_freq_1*kgrid.t_array);
   source.p(end - Nx*0.375, Ny*0.375) = source_mag*sin(2*pi*source_freq_2*kgrid.t_array);

   Is this the right way of defining two source because I keep on getting error
   'Reference to non-existent field'p'.'. Kindly advise me the right way. Your help is highly appreciated.

   Thank you

   Arslan

   Posted 10 years ago #

9. Bradley Treeby
   

   Developer
   Joined: Mar '10
   Posts: 1,643

   Hi Arslan,

   To use different time series for each source point in source.p_mask, you must define source.p as a 2D matrix. Each row of the matrix corresponds to the time series for a different source point. The ordering of the source points follows MATLAB linear matrix indexing. In your case, you have two different source points with different frequencies, so the code would be something like:

   source.p = zeros(sum(source.p_mask(:)), kgrid.Nt);
   source.p(1, :) = source_mag*sin(2*pi*source_freq_1*kgrid.t_array);
   source.p(2, :) = source_mag*sin(2*pi*source_freq_2*kgrid.t_array);

   There is some information in Section 3.4 of the k-Wave manual that might be helpful.

   Brad.

   Posted 10 years ago #

10. rarslana
    

    Member
    Joined: Jun '12
    Posts: 4

    Hi Brad!

    Thank you very much for your response its quite helpful, here I have another inquiry what if I have to add this time varying source to source mask like disc, circle etc. as initial value. DO I create source mask as

    source.p_mask = zeros(Nx, Ny);
    source.p_mask = makedisc( , , , ,);
    source.p_mask = makedisc( , , , ,);

    source.p = zeros(sum(source.p_mask(:)), kgrid.Nt);
    source.p(1, :) = source_mag*sin(2*pi*source_freq_1*kgrid.t_array);
    source.p(2, :) = source_mag*sin(2*pi*source_freq_2*kgrid.t_array);

    Kindly advise. Your cooperation is greatly appreciated.

    Regards
    Arslan

    Posted 10 years ago #

11. Bradley Treeby
    

    Developer
    Joined: Mar '10
    Posts: 1,643

    Hi Arslan,

    Again, the list of time series must follow the source mask following MATLAB's linear matrix index ordering. However, as k-Wave is currently coded, it's not completely trivial to assign different time series to different blocks of the source mask. Here's an example of how you might do it using a source mask that is defined using different regions labelled with different numbers:

    % define two different discs
    disc1 = makeDisc(Nx, Ny, Nx/4, Ny/2, 5);
    disc2 = makeDisc(Nx, Ny, 3*Nx/4, Ny/2, 4);
    
    % define two time series for each disc
    source1 = sin(2*pi*0.25e6*kgrid.t_array);
    source2 = 2*sin(2*pi*1e6*kgrid.t_array);
    
    % assign discs to source mask
    source.p_mask = zeros(Nx, Ny);
    source.p_mask(disc1 == 1) = 1;
    source.p_mask(disc2 == 1) = 2;
    
    % collapse the source mask to a 1D vector and remove the zeros
    p_mask_1D = source.p_mask(:);
    p_mask_1D(p_mask_1D == 0) = [];
    
    % use the labelled 1D source mask to assign the correct time series
    % to each source point
    source.p = zeros(sum(sum(source.p_mask == 1)), kgrid.Nt);
    source.p(p_mask_1D == 1, :) = repmat(source1, [sum(p_mask_1D == 1), 1]);
    source.p(p_mask_1D == 2, :) = repmat(source2, [sum(p_mask_1D == 2), 1]);
    
    % make the source mask binary
    source.p_mask = source.p_mask ~= 0;

    Hope that helps,

    Brad.

    Posted 10 years ago #

12. jlc4kk
    

    Member
    Joined: Jun '15
    Posts: 7

    I have a follow up question related to this last post.

    With latest k-Wave version (1.2.1), I am simulating tissue/bone heating using a transducer comprised of four concentric rings each one driven slightly off phase starting at 1 MHz, 1.001 MHz, and so on. I discretized it so there are ~50,000 active elements divided approximately evenly across the four rings. I am using kspaceFirstOrder3D with absorbion + heterogenous properties (medium.sound_speed, medium.density, and medium.alpha_coeff). In accordance with the "checkStability.m", I set the recommended time step and require about 23,000 steps total to get to steady state.

    I am running into memory constraints as my source.p is large [50,000 active elements x 23,000 num. of time steps]. This is due to copying the time series for each active element as in the example code above:

    source.p(p_mask_1D == 1, :) = repmat(source1, [sum(p_mask_1D == 1), 1]);
    source.p(p_mask_1D == 2, :) = repmat(source2, [sum(p_mask_1D == 2), 1]);

    As I only have 4 unique time series (for my 4 rings), is there a way to reduce memory here? My hope is to fit onto a GPU.

    Thank you,
    JLC

    Posted 6 years ago #

13. Bradley Treeby
    

    Developer
    Joined: Mar '10
    Posts: 1,643

    If you are using the MATLAB code, there is undocumented functionality where you can assign a non-integer source mask, where each value corresponds to a selection of grid points that belong to a particular transducer. You can then assign a time series for each transducer element (i.e., each unique positive integer in the source mask).

    As an example, open the "Monopole Point Source In A Homogeneous Propagation Medium Example" and replace the source definition with:

    % define a source points
    source.p_mask = zeros(Nx, Ny);
    source.p_mask(3*Nx/4, Ny/2:Ny/2 + 10) = 1;
    source.p_mask(Nx/4,   Ny/2:Ny/2 + 10) = 2;
    
    % define a time varying sinusoidal source
    source_freq1 = 0.25e6;   % [Hz]
    source_freq2 = 0.5e6;   % [Hz]
    source.p(1, :) = sin(2 * pi * source_freq1 * kgrid.t_array);
    source.p(2, :) = sin(2 * pi * source_freq2 * kgrid.t_array);
    
    % filter the source to remove high frequencies not supported by the grid
    source.p(1, :) = filterTimeSeries(kgrid, medium, source.p(1, :));
    source.p(2, :) = filterTimeSeries(kgrid, medium, source.p(2, :));

    Note, this functionality will be replaced with something more general in a future version (which will also work using the C++ code), so the syntax will almost certainly change.

    Hope that helps,

    Brad.

    Posted 6 years ago #

14. jlc4kk
    

    Member
    Joined: Jun '15
    Posts: 7

    Thank you, Brad. I am running the GPU code, kspaceFirstOrder3DG, from Matlab interface. Is there an estimate of when this functionality will be available?

    Posted 6 years ago #

15. Bradley Treeby
    

    Developer
    Joined: Mar '10
    Posts: 1,643

    No estimate at this stage I'm afraid. It's on the development list, just not sure if it will make it into the next release.

    Posted 6 years ago #

RSS feed for this topic

Reply

You must log in to post.