写一段超声成像的matlab代码

以下是一个简单的超声成像模拟代码，可以生成一个具有相干性的点源图像并对其进行B模式扫描：

% 超声成像模拟代码

% 设置参数
c = 1540; % 声速
fs = 40e6; % 采样频率
f0 = 5e6; % 脉冲中心频率
lambda = c/f0; % 波长
grid_size = 1e-4; % 网格大小
x = (-10:grid_size:10); % x方向网格
z = (0:grid_size:20); % z方向网格

% 生成点源图像
x0 = 0; % 点源x坐标
z0 = 5; % 点源z坐标
amp = 1; % 点源振幅
img = amp * exp(-(x-x0).^2/lambda^2) .* exp(-(z-z0).^2/lambda^2);

% 对图像进行B模式扫描
num_lines = 100; % 扫描线数量
theta = linspace(-pi/2, pi/2, num_lines); % 扫描角度
dtheta = theta(2) - theta(1); % 扫描步长
d = 10e-3; % 探头间距
width = 5e-3; % 探头宽度
t = (0:length(img)-1)/fs; % 时间轴
rf_data = zeros(length(t), num_lines); % B模式数据

for ii = 1:num_lines
    % 计算扫描线位置
    x_line = -d/2 + (ii-1)*d/(num_lines-1);
    
    % 计算扫描线响应
    for jj = 1:length(z)
        r = sqrt((x-x_line).^2 + z(jj)^2);
        t_delay = r/c;
        pulse = exp(-(t-t_delay).^2*(f0/2)^2);
        scatter_response = exp(-(x-x_line).^2/width^2);
        rf_data(jj,ii) = sum(img(:,jj).*pulse.*scatter_response);
    end
end

% 显示B模式图像
figure;
imagesc(theta*180/pi, z, abs(hilbert(rf_data)));
colormap(gray); axis image;
xlabel('扫描角度（度）');
ylabel('深度（米）');
title('B模式图像');