%===========每当打开matlab运行一下下面代码，运行一次即可，然后后注释掉。下次打开matlab再运行
% path(path,‘E:\MATLAB\kkk\Field_II\field_II_combined‘); % 加field文件夹入搜索路径
% run field_init % Initialize the Field II program system
clear;pause(0.1);
%------------------------基本参数--------------------
f0=3e6; % 中心频率[Hz]
fs=100e6; % 采样率 [Hz]
c=1540; % 声速 [m/s]
ele_width =c/f0; % 阵元宽度 x方向 使等于波长 0.51/1000
ele_height=5/1000; % 阵元高度 y方向
ele_kerf =ele_width/20; % 阵元间隙 （当width+kerf大于波长一半时，阵元间的相互作用可忽略不计。万明习（上）P228）
ele_num=64; % 阵元个数
focus=[0 0 50]/1000; % 焦点 [m]
%------------------------系统配置---------------------
set_sampling(fs); % 设置采样率
impulse_response=sin(2*pi*f0*(0:1/fs:2/f0));
impulse_response=impulse_response.*hanning(length(impulse_response))‘; % 冲击响应
excitation=sin(2*pi*f0*(0:1/fs:2/f0)); % 激励，不是单脉冲信号
emit_aperture = xdc_linear_array (ele_num, ele_width, ele_height, ele_kerf, 1, 5, focus); % 发射探头
receive_aperture= xdc_linear_array (ele_num, ele_width, ele_height, ele_kerf, 1, 5, focus); % 接收探头
xdc_impulse (emit_aperture, impulse_response); % 设置发射探头冲击响应
xdc_excitation (emit_aperture, excitation); % 设置发射探头激励
xdc_impulse (receive_aperture, impulse_response); % 设置接收探头冲击响应
%[phantom_positions, phantom_amplitudes] = cyst_phantom(100000); % 创建人工体膜，查看体膜的摆放的位置及尺寸
phantom_positions=focus;
phantom_amplitudes=1;
[RF_data, start_time]=calc_scat(emit_aperture, receive_aperture, phantom_positions, phantom_amplitudes);
plot((0:length(RF_data)-1)/fs+start_time,RF_data);
构建人工体膜，将体膜紧贴探头中心，采集RF信号
function [positions, amp] = cyst_phantom (N) % N=100000 合适
x_size = 100/1000; % Width of phantom [m]
y_size = 100/1000; % Transverse width of phantom [m]
z_size = 100/1000; % Height of phantom [m]
z_start = 0/1000; % Start of phantom surface [m];
% Create the general scatterers
x = (rand (N,1)-0.5)*x_size; % N个散射点x轴坐标均匀分布在[-25,25]
y = (rand (N,1)-0.5)*y_size; % N个散射点y轴坐标均匀分布在[-25,25]
z = rand (N,1)*z_size + z_start; % N个散射点z轴坐标均匀分布在[-25,25]+z_start
amp=rand(N,1); % N个散射点均匀分布散射强度[0,1]之间
pht_num=2; % 放置2个球
pht_r =8/1000; % 半径
dz=z_size/(pht_num+1); % 间隔
inside=zeros(N,1);
for i=1:pht_num
temp=(((x-0).^2+(y-0).^2+(z-z_start-i*dz).^2)<pht_r^2);
inside=inside|temp;
end
amp(find(inside==1))=10; % 球内的散射系数设为2
% Return the variables
positions=[x y z];
% 只绘制球内的散点
figure;
for i=1:N
if inside(i)
plot3(x(i),y(i),z(i),‘o‘);hold on;
end
end
axis([-50/1000 50/1000 -50/1000 50/1000 0 100/1000]);
grid on;xlabel(‘x‘);ylabel(‘y‘);zlabel(‘z‘);
end