k-Wave User Forum » Topic: How to obtain harmonic components
http://www.k-wave.org/forum/topic/how-to-obtain-harmonic-components
Support for the k-Wave MATLAB toolboxen-USMon, 05 Aug 2024 09:45:11 +0000http://bbpress.org/?v=1.0.2<![CDATA[Search]]>q
http://www.k-wave.org/forum/search.php
jackYANG on "How to obtain harmonic components"
http://www.k-wave.org/forum/topic/how-to-obtain-harmonic-components#post-7670
Sun, 28 Jun 2020 02:42:26 +0000jackYANG7670@http://www.k-wave.org/forum/<p>Hi Brad,<br />
Thanks for your quick reply. I saw it early this morning.I want to make it clear that the first problem is actually what I saw in example_us_beam_patterns. M from the Kwave example, and the details are as follows:<br />
% reshape the sensor data to its original position so that it can be<br />
% indexed as sensor_data(x, j, t)<br />
sensor_data = reshape(sensor_data, [Nx, Nj, kgrid.Nt]);</p>
<p>% compute the amplitude spectrum<br />
[freq, amp_spect] = spect(sensor_data, 1/kgrid.dt, 'Dim', 3);</p>
<p>% compute the index at which the source frequency and its harmonics occur<br />
[f1_value, f1_index] = findClosest(freq, tone_burst_freq);<br />
[f2_value, f2_index] = findClosest(freq, 2 * tone_burst_freq);</p>
<p>% extract the amplitude at the source frequency and store<br />
beam_pattern_f1 = amp_spect(:, :, f1_index);</p>
<p>% extract the amplitude at the second harmonic and store<br />
beam_pattern_f2 = amp_spect(:, :, f2_index); </p>
<p>% extract the integral of the total amplitude spectrum<br />
beam_pattern_total = sum(amp_spect, 3);</p>
<p>% plot the beam patterns<br />
figure;<br />
imagesc(j_vec * 1e3, (kgrid.x_vec - min(kgrid.x_vec(:))) * 1e3, beam_pattern_f1 * 1e-3);<br />
xlabel([j_label '-position [mm]']);<br />
ylabel('x-position [mm]');<br />
title('Beam Pattern At Source Fundamental');<br />
colormap(jet(256));<br />
c = colorbar;<br />
ylabel(c, 'Pressure [kPa]');<br />
axis image;</p>
<p>figure;<br />
imagesc(j_vec * 1e3, (kgrid.x_vec - min(kgrid.x_vec(:))) * 1e3, beam_pattern_f2 * 1e-3);<br />
xlabel([j_label '-position [mm]']);<br />
ylabel('x-position [mm]');<br />
title('Beam Pattern At Second Harmonic');<br />
colormap(jet(256));<br />
c = colorbar;<br />
ylabel(c, 'Pressure [kPa]');<br />
axis image;</p>
<p>figure;<br />
imagesc(j_vec * 1e3, (kgrid.x_vec - min(kgrid.x_vec(:))) * 1e3, beam_pattern_total * 1e-6);<br />
xlabel([j_label '-position [mm]']);<br />
ylabel('x-position [mm]');<br />
title('Total Beam Pattern Using Integral Of Recorded Pressure');<br />
colormap(jet(256));<br />
c = colorbar;<br />
ylabel(c, 'Pressure [MPa]');<br />
axis image;</p>
<p>% =========================================================================<br />
% PLOT DIRECTIVITY PATTERN AT FOCUS<br />
% =========================================================================</p>
<p>% compute the directivity at each of the harmonics<br />
directivity_f1 = squeeze(beam_pattern_f1(round(transducer.focus_distance/dx), :));<br />
directivity_f2 = squeeze(beam_pattern_f2(round(transducer.focus_distance/dx), :));</p>
<p>% normalise<br />
directivity_f1 = directivity_f1 ./ max(directivity_f1(:));<br />
directivity_f2 = directivity_f2 ./ max(directivity_f2(:));</p>
<p>% compute relative angles from transducer<br />
if strcmp(MASK_PLANE, 'xy')<br />
horz_axis = ((1:Ny) - Ny/2) * dy;<br />
else<br />
horz_axis = ((1:Nz) - Nz/2) * dz;<br />
end<br />
angles = 180 * atan2(horz_axis, transducer.focus_distance) / pi;</p>
<p>% plot the directivity<br />
figure;<br />
plot(angles, directivity_f1, 'k-', angles, directivity_f2, 'k--');<br />
axis tight;<br />
xlabel('Angle [deg]');<br />
ylabel('Normalised Amplitude');<br />
legend('Fundamental', 'Second Harmonic', 'Location', 'NorthWest'); </p>
<p>For the second question, I don't know how to upload photos, please forgive me.<br />
Is my with phased array (in water) emission is focused on the sound waves, respectively by fat and muscle, focus of sound waves in the waters, result in focus position I saw there was a shock (axial and radial pressure curve, the focus position there will be a small dent or bulging xiao feng), I am using p_max or p_min records of axial pressure curve, has this kind of phenomenon, if I use p_rms is very smooth.<br />
I don't quite understand,Hope to get your answer.<br />
jackYANG
</p>Bradley Treeby on "How to obtain harmonic components"
http://www.k-wave.org/forum/topic/how-to-obtain-harmonic-components#post-7669
Sat, 27 Jun 2020 20:20:53 +0000Bradley Treeby7669@http://www.k-wave.org/forum/<p>Hi Jack,</p>
<p>Unless I'm missing something, your code (1) above does not compute the root-mean-square of the pressure field. In general though, the <code>p_rms</code> is probably not the most useful measure if you are simulating nonlinear fields. Another unrelated point, I'm assuming you have reshaped <code>sensor_data.p</code>. If not, the time dimension when using a binary or Cartesian sensor mask is dim 2 not dim 3.</p>
<p>I don't understand your second question I'm afraid.</p>
<p>Brad.
</p>jackYANG on "How to obtain harmonic components"
http://www.k-wave.org/forum/topic/how-to-obtain-harmonic-components#post-7664
Sat, 27 Jun 2020 14:16:16 +0000jackYANG7664@http://www.k-wave.org/forum/<p>Hi brad,<br />
I'm excited to see you answering our learner's questions.I would like to ask what is the difference between the two kinds of calculation of RMS sound pressure, I clearly see in my simulation results that these two kinds of calculation results are very different,I would like to ask what is the difference between the two kinds of calculation of RMS sound pressure, I clearly see in my simulation results that these two kinds of calculation results are very different.<br />
1<br />
[freq, amp_spect] = spect(sensor_data.p, 1/kgrid.dt, 'Dim', 3);<br />
[f11_value, f11_index] = findClosest(freq, tone_burst_freq);<br />
[f12_value, f12_index] = findClosest(freq, 2 * tone_burst_freq);<br />
beam_pattern_total = sum(amp_spect, 3);<br />
2<br />
sensor.record={'p_rms'};</p>
<p>I would like to ask you another question, why I simulate the acoustic axis (radial) sound pressure curve with P_max and P_min after the phased array focused ultrasound passes through two layers of tissue (fat and muscle), and the focus will be oscillated, as if the waveform is distorted
</p>jackYANG on "How to obtain harmonic components"
http://www.k-wave.org/forum/topic/how-to-obtain-harmonic-components#post-7610
Sat, 13 Jun 2020 14:00:44 +0000jackYANG7610@http://www.k-wave.org/forum/<p>Ps:f0=1MHz,fmax=5MHz,I am sure that there is a strong nonlinear change at the focal point and that the spectrum analysis at the focal point yields a normal distribution (1-5mhz)
</p>jackYANG on "How to obtain harmonic components"
http://www.k-wave.org/forum/topic/how-to-obtain-harmonic-components#post-7609
Sat, 13 Jun 2020 13:57:15 +0000jackYANG7609@http://www.k-wave.org/forum/<p>Hi Brad,<br />
This is my code. I hope you can take some time to help me analyze it. Thank you very much for your help and kindness in learning Kwave.<br />
clc;<br />
clear all;<br />
close all;<br />
DATA_CAST = 'single'; % set to 'single' or 'gpuArray-single' to speed up computations<br />
MASK_PLANE = 'xz'; % set to 'xy' or 'xz' to generate the beam pattern in different planes<br />
USE_STATISTICS = true; % set to true to compute the rms or peak beam patterns, set to false to compute the harmonic beam patterns<br />
% =========================================================================<br />
% SIMULATION<br />
% =========================================================================</p>
<p>% 创建完美匹配层(PML)<br />
PML_X_SIZE = 10; % [grid points]<br />
PML_Y_SIZE = 10; % [grid points]<br />
PML_Z_SIZE = 10; % [grid points]</p>
<p>Nx = 430; % number of grid points in the x (row) direction<br />
Ny =430; % number of grid points in the y (column) direction<br />
Nz =430;<br />
% x_size=72e-3;<br />
dx = 0.15e-3; % grid point spacing in the x direction [m]<br />
dy = dx; % grid point spacing in the y direction [m]<br />
dz = dx;<br />
kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz);<br />
% 定义传播介质的属性</p>
<p>medium.sound_speed = 1500*ones(Nx,Ny,Nz) ; % [m/s] 2nd<br />
medium.alpha_coeff = 0.0022*ones(Nx,Ny,Nz); % [dB/(MHz^y cm)] 2nd<br />
medium.density = 1000*ones(Nx,Ny,Nz); % [kg/m^3] 2nd</p>
<p>% medium.sound_speed(:,:,60:70) = 1615 ; % [m/s] 1st 101<br />
% medium.density(:,:,60:70) = 1090; % [kg/m^3] 1st<br />
% medium.alpha_coeff(:,:,60:70)=1.8 ; % [dB/(MHz^y cm)] 1st</p>
<p>medium.sound_speed(:,:,100:300) = 1479 ; % [m/s] 1st 101<br />
medium.density(:,:,100:300) = 928; % [kg/m^3] 1st<br />
medium.alpha_coeff(:,:,100:300)=0.6 ; % [dB/(MHz^y cm)] 1st<br />
%<br />
medium.sound_speed(:,:,301:end) = 1588 ; % [m/s] 1st 101<br />
medium.density(:,:,301:end) = 1130; % [kg/m^3] 1st<br />
medium.alpha_coeff(:,:,301:end)=0.7; % [dB/(MHz^y cm)] 1st<br />
%<br />
medium.alpha_power = 2;<br />
medium.BonA=zeros(Nx,Ny,Nz);<br />
medium.BonA = 5*ones(Nx,Ny,Nz);<br />
medium.BonA(:,:,100:300) =10;<br />
medium.BonA(:,:,301:end)=7.8;<br />
% medium.BonA=18*ones(Nx,Ny,Nz);<br />
% kgrid.t_array = 0:2e-8:2.352e-5;<br />
% dt=2e-8;<br />
% Nt=1176;<br />
% kgrid.t_array=0:dt:Nt;<br />
kgrid.makeTime(medium.sound_speed);</p>
<p>%% creat element<br />
diameter = 60*1e-3;<br />
radius = 60*1e-3;<br />
center_pos = [1,1,0.05]*(dx*Nx)/2;<br />
num_per_loop = [4,8,12,16,24];<br />
number_element = sum(num_per_loop);<br />
bowl_pos = makeCartLoopBowl(diameter,radius, center_pos,num_per_loop, true);</p>
<p>%% mask setting<br />
source.p_mask = zeros(Nx,Ny,Nz);<br />
bowl_pos_n = round(bowl_pos/dx);<br />
for j = 1:number_element<br />
source.p_mask(bowl_pos_n(1,j),bowl_pos_n(2,j),bowl_pos_n(3,j))=1;<br />
end</p>
<p>%% create the input signal using toneBurst<br />
source_strength = 1e6; % [Pa]<br />
tone_burst_freq = 1e6 ; % [Hz]<br />
tone_burst_cycles = 5;<br />
input_signal = source_strength*toneBurst(1/kgrid.dt, tone_burst_freq, tone_burst_cycles,'SignalLength',kgrid.Nt);<br />
source.p = input_signal;<br />
% source.p_mode='additive';<br />
source.p_mode='dirichlet';<br />
%% caltuate phase_control delay time<br />
sound_speed =1500;<br />
focus_position = [0 ,0, 0.0294]; %笛卡尔坐标表示<br />
input_signal_mat = focus(kgrid, source.p, source.p_mask, focus_position, sound_speed);<br />
row = find(source.p_mask==1);<br />
source.p = input_signal_mat;</p>
<p>%% creat area elemrnt<br />
bowl_pos = round(bowl_pos'/dx);<br />
radius = 400;<br />
diameter = 41;<br />
focus_pos = [216, 216, 412];<br />
grid_size = [Nx,Ny,Nz];<br />
[bowl_pos,bowls_labelled] = makeMultiBowl(grid_size, bowl_pos, radius, diameter, focus_pos, 'Plot', true);<br />
%source.p_mask = zeros(Nx, Ny,Nz);<br />
%[grid_data, order_index, reorder_index] = cart2grid(kgrid, sensor.mask);<br />
source.p_mask=bowl_pos;</p>
<p>%% apply the signal in the area element<br />
num = find(source.p_mask==1);<br />
input_signal_mat = input_signal_mat(:,1:kgrid.Nt);<br />
input_signal_new = [];<br />
input_signal_spread = zeros(length(num),kgrid.Nt);<br />
for j = 1:number_element<br />
ind = find(bowls_labelled==j);<br />
index_element = length(ind);<br />
index_row = [];<br />
for i = 1:index_element<br />
index_label = find(abs(row - ind(i))<=2);<br />
index_row = [index_row;index_label];<br />
end</p>
<p> index_row = index_row(1);<br />
input_signal_new = [input_signal_new;input_signal_mat(index_row,:)];<br />
end<br />
for j = 1:number_element<br />
ind = find(bowls_labelled==j);<br />
index_element = length(ind);<br />
for i = 1:index_element<br />
cal = find(abs(num-ind(i))<1);<br />
input_signal_spread(cal(1),:) = input_signal_new(j,:);<br />
end<br />
end<br />
source.p = input_signal_spread;<br />
% source.p_mode='additive';<br />
source.p_mode='dirichlet';<br />
sensor.mask= zeros(Nx, Ny,Nz);</p>
<p>switch MASK_PLANE<br />
case 'xy'</p>
<p> % 定义掩膜<br />
sensor.mask(:, :, Nz/2+1) = 1;</p>
<p> % 存储y轴属性<br />
Nj = Ny;<br />
j_vec = kgrid.y_vec;<br />
j_label = 'y';</p>
<p> case 'xz'</p>
<p> % 定义掩膜<br />
sensor.mask(:, Ny/2+1, :) = 1;</p>
<p> % 存储z轴属性<br />
Nj = Nz;<br />
j_vec = kgrid.z_vec;<br />
j_label = 'z';</p>
<p>end<br />
input_args = {'PMLInside', false, 'PlotPML', false, 'PMLSize', [PML_X_SIZE, PML_Y_SIZE, PML_Z_SIZE],'PlotLayout',false, 'Smooth',[true,true,true], 'PlotPML', false,'PlotScale',[-source_strength,source_strength]};<br />
% run the simulation<br />
sensor.record={'p_max','p','p_rms','p_final','p_min'};<br />
sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});<br />
sensor_data.p= reshape(sensor_data.p, [Nx, Nz,kgrid.Nt]);<br />
p1 = zeros(1,kgrid.Nt);<br />
p1 = reshape(sensor_data.p(216,412,:),[1,kgrid.Nt]);<br />
spect(p1,1/kgrid.dt,'Plot',true)</p>
<p>if USE_STATISTICS</p>
<p> sensor_data.p_rms = reshape(sensor_data.p_rms, [Nx, Nz]);<br />
sensor_data.p_max = reshape(sensor_data.p_max, [Nx, Nz]);<br />
sensor_data.p_min = reshape(sensor_data.p_min, [Nx, Nz]);</p>
<p> figure(1);<br />
imagesc(j_vec * 1e3, (kgrid.x_vec - min(kgrid.x_vec(:))) * 1e3, sensor_data.p_max * 1e-6);<br />
xlabel([j_label '-position [mm]']);<br />
ylabel('x-position [mm]');<br />
title('Total Beam Pattern Using Maximum Of Recorded Pressure');<br />
colormap(jet(256));<br />
c = colorbar;<br />
ylabel(c, 'Pressure [MPa]');<br />
axis image</p>
<p> figure(2);<br />
imagesc(j_vec * 1e3, (kgrid.x_vec - min(kgrid.x_vec(:))) * 1e3, sensor_data.p_rms * 1e-6);<br />
xlabel([j_label '-position [mm]']);<br />
ylabel('x-position [mm]');<br />
title('Total Beam Pattern Using RMS Of Recorded Pressure');<br />
colormap(jet(256));<br />
c = colorbar;<br />
ylabel(c, 'Pressure [MPa]');<br />
axis image;</p>
<p> figure(3);<br />
imagesc(j_vec * 1e3, (kgrid.x_vec - min(kgrid.x_vec(:))) * 1e3, sensor_data.p_min * 1e-6);<br />
xlabel([j_label '-position [mm]']);<br />
ylabel('x-position [mm]');<br />
title('Total Beam Pattern Using p_min Of Recorded Pressure');<br />
colormap(jet(256));<br />
c = colorbar;<br />
ylabel(c, 'Pressure [MPa]');<br />
axis image;</p>
<p> return<br />
end
</p>Bradley Treeby on "How to obtain harmonic components"
http://www.k-wave.org/forum/topic/how-to-obtain-harmonic-components#post-7608
Sat, 13 Jun 2020 13:28:15 +0000Bradley Treeby7608@http://www.k-wave.org/forum/<p>Hi jackYANG,</p>
<p>It's impossible to say without studying the problem in some detail, which unfortunately I don't have capacity to do. The only things I can suggest are that you (1) setup a simple scenario and compare with analytical solutions to make sure you are using k-Wave correctly, and (2) run a proper convergence test (e.g., decreasing the grid spacing and CFL) to make sure you can trust your numerical results for the more complicated case.</p>
<p>Brad.
</p>jackYANG on "How to obtain harmonic components"
http://www.k-wave.org/forum/topic/how-to-obtain-harmonic-components#post-7607
Sat, 13 Jun 2020 13:24:48 +0000jackYANG7607@http://www.k-wave.org/forum/<p>Hi Brad,<br />
I hope you can see the question.I recently reduced the DX by a factor of 0.15E-3m, equivalent to one-tenth of the wavelength. I used FCous function to conduct the in-water focus under linear conditions, and the focus shift position was 0.3mm. When I added tissue, the focus shift position was 0.8mm,do you find it credible?
</p>jackYANG on "How to obtain harmonic components"
http://www.k-wave.org/forum/topic/how-to-obtain-harmonic-components#post-7577
Sun, 07 Jun 2020 03:58:21 +0000jackYANG7577@http://www.k-wave.org/forum/<p>Hi Brad,<br />
First of all, thank you for your answer.In fact,I am trying to exactly replicate a published paper,but the approach is slightly different,in the published paper is FDTD,but I wanted to do it with the KWAVE (PSTD),then observe the effect of some parameters on the focus shift and sound field method,But I did change a lot of the parameters and didn't get a deviation of around 3mm (actually only around 0.3mm).Other sound field patterns seem to be consistent.I'm now running the reduced dx case to allow more harmonics to pass through, hopefully to see what I want to see, by the way, do you have any questions or Suggestions about the code that I sent you in the program.Finally, thank you for your patience and for your guidance.</p>
<p> jackYANG
</p>Bradley Treeby on "How to obtain harmonic components"
http://www.k-wave.org/forum/topic/how-to-obtain-harmonic-components#post-7566
Sat, 06 Jun 2020 18:08:25 +0000Bradley Treeby7566@http://www.k-wave.org/forum/<p>HI jackYANG,</p>
<p>Any shifts in the focus will depend on many factors, including the exact geometry, sound speed, level of nonlinearity, etc. You can definitely use k-Wave to simulate the shift in focus due to nonlinearity, although the shift is very small if the nonlinearity is weak (your reported grid parameters only let you simulate the second harmonic). Are you trying to exactly replicate a published paper, or just observe these shifts in general? </p>
<p>Brad.
</p>jackYANG on "How to obtain harmonic components"
http://www.k-wave.org/forum/topic/how-to-obtain-harmonic-components#post-7531
Wed, 27 May 2020 15:10:27 +0000jackYANG7531@http://www.k-wave.org/forum/<p>hello Brad and the others,<br />
I really hope someon e can help me understand and solve this problem,I've been working on this problem for a few days,my biological tissue(fat and muscle) parameters are as follows:<br />
medium.sound_speed = 1500*ones(Nx,Ny,Nz) ; % [m/s] 2nd<br />
medium.alpha_coeff = 0.0022*ones(Nx,Ny,Nz); % [dB/(MHz^y cm)] 2nd<br />
medium.density = 1000*ones(Nx,Ny,Nz); % [kg/m^3] 2nd</p>
<p>medium.sound_speed(:,:,80:150) = 1478 ; % [m/s] 1st 101<br />
medium.density(:,:,80:150) = 925; % [kg/m^3] 1st<br />
medium.alpha_coeff(:,:,80:150)=0.6 ; % [dB/(MHz^y cm)] 1st</p>
<p>end=240,dx=0.3e-3;<br />
medium.sound_speed(:,:,151:end) = 1588 ; % [m/s] 1st<br />
medium.density(:,:,151:end) = 1090; % [kg/m^3] 1st<br />
medium.alpha_coeff(:,:,151:end)=0.7 ; % [dB/(MHz^y cm)] 1st</p>
<p>medium.alpha_power = 2;<br />
medium.BonA=6;<br />
My transducer is a concave sphere with multiple elements(64),R=D=60mm,r=d=21*dx;<br />
When I use the Fcous implementation to focus through these biological tissues,<br />
the actual focus location overlaps with the default fcous_pos.<br />
When I change the speed of sound of fat from 1478 to 3000, t<br />
he focus position will have a big deviation of 1.3cm.I learned from relevant literature that through this fcous focusing method, the focus is generally shifted by about 0.3cm.</p>
<p>PS:dt=60ns,Nt=1500;</p>
<p>Another question,when I multiply media.alpha_coeff for each layer (except for the first layer) by 2 or 4,medium.alpha_power no change,according to the literature, the focus is also slightly shifted forward,but I haven't seen such a change in the kwave simulation so far,anybody know why?
</p>jackYANG on "How to obtain harmonic components"
http://www.k-wave.org/forum/topic/how-to-obtain-harmonic-components#post-7516
Thu, 21 May 2020 12:37:10 +0000jackYANG7516@http://www.k-wave.org/forum/<p>Hi Brad,<br />
Thank you for your response.I will have a try.Now another problem began to haunt me,I used the foucs function in Kwave to focus nonlinearly on multiple layers of tissue,theoretically, due to refraction, absorption, and nonlinearity, the desired focus is shifted,however, in my simulation, a strange phenomenon occurs, that is, the actual focus and the expected focus are basically in the same position,I changed a lot of parameters, but none of them allowed me to see the shift in focus,the concave spherical multielement phased array I used, f0=1MHz,fmax= 2.5mhz.I'm looking forward to your quick answer. It's so important to me.<br />
sincerrly<br />
JackYANG
</p>Bradley Treeby on "How to obtain harmonic components"
http://www.k-wave.org/forum/topic/how-to-obtain-harmonic-components#post-7513
Thu, 21 May 2020 11:11:21 +0000Bradley Treeby7513@http://www.k-wave.org/forum/<p>Hi Jack,</p>
<p>Unfortunately, if you want to simulate high frequencies in a large domain, you need a large grid. Currently, there is no way in k-Wave to use a non-uniform grid. If your transducer and medium are axisymmetric, using the axisymmetric code will significantly reduce the computational load.</p>
<p>Brad.
</p>jackYANG on "How to obtain harmonic components"
http://www.k-wave.org/forum/topic/how-to-obtain-harmonic-components#post-7504
Mon, 18 May 2020 03:12:07 +0000jackYANG7504@http://www.k-wave.org/forum/<p>Hi,Brad,<br />
Sorry to bother you again,the simulation I'm doing right now has a wide range of domains,I want to capture more harmonic components that propagate through the grid,but this will cause my grid to be very large (2048*2048*2048),this would be an impossible task on my computer,I have tried using gpu for simulation before,but it failed,can you tell me how to do it? Can you give me a simple actual example?In addition,I wonder in what way can I define a non-uniform grid to capture more harmonic component simulations,I hope to get your reply.
</p>