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
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
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
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
Hi Chao,
It shouldn't do, if the absorption's turned off. Can you post your code?
Thanks,
Ben
Hi Dr. Cox,
Just ignore my last post. The amplitude does not decay. It's my mistake. Sorry about that.
Chao
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
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.
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
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.
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
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.
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
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.
Thank you, Brad. I am running the GPU code, kspaceFirstOrder3DG, from Matlab interface. Is there an estimate of when this functionality will be available?
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.
You must log in to post.