http://www.k-wave.org/forum/topic/how-to-obtain-harmonic-components Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 00:19:04 +0000 http://bbpress.org/?v=1.0.2 q http://www.k-wave.org/forum/search.php http://www.k-wave.org/forum/topic/how-to-obtain-harmonic-components#post-7670 Sun, 28 Jun 2020 02:42:26 +0000 jackYANG 7670@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> http://www.k-wave.org/forum/topic/how-to-obtain-harmonic-components#post-7669 Sat, 27 Jun 2020 20:20:53 +0000 Bradley Treeby 7669@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> http://www.k-wave.org/forum/topic/how-to-obtain-harmonic-components#post-7664 Sat, 27 Jun 2020 14:16:16 +0000 jackYANG 7664@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> http://www.k-wave.org/forum/topic/how-to-obtain-harmonic-components#post-7610 Sat, 13 Jun 2020 14:00:44 +0000 jackYANG 7610@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> http://www.k-wave.org/forum/topic/how-to-obtain-harmonic-components#post-7609 Sat, 13 Jun 2020 13:57:15 +0000 jackYANG 7609@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))&lt;=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))&lt;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> http://www.k-wave.org/forum/topic/how-to-obtain-harmonic-components#post-7608 Sat, 13 Jun 2020 13:28:15 +0000 Bradley Treeby 7608@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> http://www.k-wave.org/forum/topic/how-to-obtain-harmonic-components#post-7607 Sat, 13 Jun 2020 13:24:48 +0000 jackYANG 7607@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> http://www.k-wave.org/forum/topic/how-to-obtain-harmonic-components#post-7577 Sun, 07 Jun 2020 03:58:21 +0000 jackYANG 7577@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> http://www.k-wave.org/forum/topic/how-to-obtain-harmonic-components#post-7566 Sat, 06 Jun 2020 18:08:25 +0000 Bradley Treeby 7566@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> http://www.k-wave.org/forum/topic/how-to-obtain-harmonic-components#post-7531 Wed, 27 May 2020 15:10:27 +0000 jackYANG 7531@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> http://www.k-wave.org/forum/topic/how-to-obtain-harmonic-components#post-7516 Thu, 21 May 2020 12:37:10 +0000 jackYANG 7516@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> http://www.k-wave.org/forum/topic/how-to-obtain-harmonic-components#post-7513 Thu, 21 May 2020 11:11:21 +0000 Bradley Treeby 7513@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> http://www.k-wave.org/forum/topic/how-to-obtain-harmonic-components#post-7504 Mon, 18 May 2020 03:12:07 +0000 jackYANG 7504@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>