Hi I'm working with k-wave
First I did a Delay & sum image reconstruction with Phased Array in 2D simulation
so i want to adopt a nonlinear effects with this Phased array
make a half of transducers in array was exictied in different frequency and
make a nonlinear effect at focal points but
when i did in nonlinear examples it doesn't work so I wonder what is the problem
here it is my code
clear all;
% =========================================================================
% DEFINE SIMULATION PROPERTIES
% =========================================================================
% define the properties used in the simulation
p0 = 7*10^6; % source pressure [Pa]
c0 = 1500; % sound speed [m/s]
rho0 = 1000; % density [kg/m^3]
alpha_0 = 0.25; % absorption coefficient [dB/(MHz^2 cm)]
sigma = 1; % shock parameter
source_freq = 1e6; % frequency [Hz]
points_per_wavelength = 50; % number of grid points per wavelength at f0
wavelength_separation = 15; % separation between the source and detector
pml_size = 64; % PML size
CFL = 0.25; %
% compute corresponding grid spacing
dx = c0/(points_per_wavelength*source_freq); % [m]
dy = c0/(points_per_wavelength*source_freq);
% compute corresponding grid size
Nx = wavelength_separation*points_per_wavelength + 20;
Ny = wavelength_separation*points_per_wavelength + 20;
% =========================================================================
% RUN SIMULATION
% =========================================================================
% create the computational grid
kgrid = makeGrid(Nx, dx, Ny, dy );
% assign the properties of the propagation medium
medium.sound_speed = ones(Nx,Ny)*c0;
medium.density = ones(Nx,Ny)*rho0;
medium.alpha_power = 2;
medium.alpha_coeff = alpha_0;
medium.alpha_mode = 'no_dispersion';
% extract the maximum frequency supported by the grid
f_max = min(min(medium.sound_speed))/(2*dx); % [Hz]
Num_T = 128;
% define a single source element
source.p_mask = zeros(Nx,Ny);
source.p_mask(Nx,Ny/2-Num_T+2:+2:Ny/2+Num_T) = 1;
% define a single sensor position an integer number of wavelengths away
sensor.mask = zeros(Nx, Ny);
sensor.mask(Nx/2,Ny/2)=1;
x_px = wavelength_separation*points_per_wavelength;
x = x_px*dx;
% compute the nonlinearity coefficient required to give the correct shock
% parameter
mach_num = p0/(rho0*c0.^2);
k = 2*pi*source_freq/c0;
BonA = 2*(sigma/(mach_num*k*x) - 1);
medium.BonA = BonA;
display_mask = source.p_mask;
% set the simulation options
input_args = {'PMLInside', false, 'PMLSize', [pml_size,pml_size], 'PMLAlpha', 1.5,'DisplayMask', display_mask...
,'RecordMovie', true, 'MovieType', 'image', 'MovieName', 'Non_linear'};
% compute points per temporal period
points_per_period = round(points_per_wavelength / CFL);
% compute corresponding time spacing
dt = 1/(points_per_period*source_freq);
focus_X=Nx/2*dx;
focus_Y=0;
for i=1:1:Num_T/2
tone_burst_offset(1,i+Num_T/2)=(sqrt(focus_X^2+focus_Y^2)-sqrt(((2*i-1)*dy-focus_Y)^2+focus_X^2))/1500/dt ;
end
for i=-1:-1:-Num_T/2
tone_burst_offset(1,i+Num_T/2+1)=(sqrt(focus_X^2+focus_Y^2)-sqrt(((2*i+1)*dy-focus_Y)^2+focus_X^2))/1500/dt ;
end
Min=min(tone_burst_offset(1,:));
for i=1:1:Num_T
tone_burst_offset(1,i)= tone_burst_offset(1,i)-Min;
end
tone_burst_freq = source_freq*1; % [Hz]
tone_burst_cycles = 10;
source.p = toneBurst(1/dt, tone_burst_freq, tone_burst_cycles, 'SignalOffset', tone_burst_offset)
source.p = (p0./(1500*1000)).*source.p;
source_lf = toneBurst(1/dt, tone_burst_freq+10^5, tone_burst_cycles, 'SignalOffset', tone_burst_offset);
source_lf = (p0./(1500*1000)).*source.p*2;
for a=1:Num_T/2
for b=1:length(source.p)
source.p(a,b) = source_lf(a,b);
end
end
% create the time array using an integer number of points per period
t_end = Nx*dx/1500;
kgrid.t_array = 0:dt:t_end;
% create the source term, offset by dt/2 so there is a point that lands on the axis
% source.p = p0*sin(2*pi*source_freq*(kgrid.t_array));
% run the simulation
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});
plot(sensor_data)
save 'data_set'
[f, sensor_data_as] = spect(sensor_data(1, :), 1/dt);
figure;
[f_sc, scale, prefix] = scaleSI(max(f(:)));
plot(scale*f, sensor_data_as, 'k-');
ylabel('Amplitude Spectrum');
xlabel(['Frequency [' prefix 'Hz]']);
legend('Original Time Series', 'Filtered Time Series');
set(gca, 'XLim', [0 1]);