不是,是超声相控阵成像Nt = kgrid.Nt;
dt = kgrid.dt;
e_tv_theta = zeros(Nt, N_elements, N_waves);
for i = 1:N_elements
for j = 1:N_waves
e_tv_theta(:,i,j) = filtered_element_data(:, i, j);
end
end
% PARAMETERS
fs = 1/kgrid.dt;
% f = tone_burst_freq;
fLow = 3000000; % [Hz]3.5MHz
fHigh = 7000000; % [Hz]6.5MHz
% c1 = c / 2;
% pitch = transducer.element_width + transducer.element_spacing; % [grid points]
% v = round(prb_x0./dx) + Nx/2;
v = 1:N_elements;
% v = u; % Receive
[nt,nv,nw] = size(e_tv_theta); % axis
zStart = 0; % longitudinal
zEnd = grid_depth; % [m]
dz = dx;
vStep = pitch; % [m]u步长(m)
% --------------------
% Spectral coordinate
% --------------------
nFFTt = 2^nextpow2(nt);
nFFTv = 2^nextpow2(nv);
% nFFTv = 2^nextpow2(nv);
% f1 = ((0:(nt-1)) - floor(nt/2)+0.5)'/((nt-1)*dt); % [Hz]中心化的频率轴
% omega = ifftshift(f1)*2*pi;
f1 = ((0:(nFFTt-1)) - floor(nFFTt/2)+0.5)'/nFFTt/dt; % [Hz]中心化的频率轴
% f1 = (0:nFFTt/2)'/nFFTt/dt;
omega = f1*2*pi;
omega = ifftshift(omega);
omegaBand = (omega >=-2*pi*fLow)& (omega <= 2*pi*fHigh);
% omegaBand = (omega >=-2*pi*fLow)& (omega <= 2*pi*fHigh); % 角频率边界限制
% omega=omega(omegaBand);
% k= omega(omegaBand)/c;
% kvs = vStep; % 波数域步长
% kv = ((0:(nv-1)) - floor(nv/2)+0.5)/((nv-1)*kvs); % 波数域变量列表
% kv = ifftshift(kv);
% kv = ku;
kvs = (2*pi)/vStep; % 波数域步长
kv = ((0:(nFFTv-1)) - floor(nFFTv/2)+0.5)*(kvs/nFFTv); % 波数域变量列表
kv = ifftshift(kv);
% E_wkv_theta = fft2(e_tv_theta,nFFTt,nFFTv);
% signal=e_tv_theta(:,32,71);t = (0:length(signal)-1) * kgrid.dt;
% figure;plot(t,signal);
% signal=E_wkv_theta(:,32,3);figure;plot(f1,fftshift(abs(signal)));
% E_wkv_theta = fft(fft(e_tv_theta,nFFTt,1),nFFTv,2);
% E_wkv_theta = E_wkv_theta(omegaBand,:,:); % e(t,v,theta)→E(ω,ku,theta)
% E_wuv = fft(e_tvw,nFFTt,1);
% --------------------
% Stolt
% --------------------
% KX = KU + KV; % kx = ku + kv
% KZ = sqrt(k.^2 - KU.^2) + sqrt(k.^2 - KV.^2);
nkx = round((N_elements-1)*pitch/dx);
% nkx = Nx;
kxs = (2*pi)/dx;
% nkx = N_elements;
% kxs = (2*pi)/pitch;
kx = ((0:(nkx-1)) - floor(nkx/2)+0.5)*(kxs/nkx);
kx = ifftshift(kx);
nkz = Nz;
kzs = kxs;
kz = ((0:(nkz-1)) - floor(nkz/2)+0.5)*(kzs/nkz);
% kz = (0:nkz/2)*(kzs/nkz);
kz = ifftshift(kz);
[KZ, KX] = ndgrid(kz, kx);
G_kz_kx = zeros(nkz, nkx);
% c1 = c0/2;
[Omega_ori, Kv_ori] = ndgrid(omega(omegaBand), kv);
for iw = 1:N_waves
e_tv = squeeze(e_tv_theta(:, :, iw));
E_wkv = fft2(e_tv,nFFTt,nFFTv);
E_wkv = E_wkv(omegaBand,:);
% temp=delay_values(iw,:)-delay(iw);
temp=0.5*grid_width*sin(angles(iw))/c0;
E_wkv_com = E_wkv.*exp(1i.*Omega_ori.*temp);
% E_wkv_com = E_wkv.*exp(1i.*Omega_ori.*delay(3));
% .*exp(i.*Omega_ori*delay(iw));
% E_wkv = squeeze(E_wkv_theta(:, :, iw));
% signal=E_wkv(:,32);figure;plot(f1,fftshift(abs(signal)));
% E(ω,kv) 每个theta对应ω, kv 网格的复数据
% t = ((prb_x0(N_elements) + grid_width/2).*sin(angles(iw)))/c0;
% A_w = exp(1i.*Omega_ori.* t);
% E_wkv = E_wkv.*A_w;
%
kx_target = KX;
kz_target = KZ;
k_target = (kx_target.^2 + kz_target.^2)./(2.*kx_target.*sin(angles(iw)) + 2.* kz_target.*cos(angles(iw))+ eps);
omega_target = c0 .* abs(k_target);
invalid_mask = isnan(omega_target) | ~isreal(omega_target);
omega_target(invalid_mask) = 0;
kv_target = kx_target - k_target .* sin(angles(iw));
valid_mask = (k_target > 0) & (k_target.^2 >= kv_target.^2) & ...
(omega_target >= min(omega(omegaBand))) & ...
(omega_target <= max(omega(omegaBand))) & ...
(abs(kv_target) <= max(abs(kv)));
% 插值 E(ω,kv)→E(kx,kz)
F_interp = griddata(...
Omega_ori, Kv_ori, E_wkv_com, ...
omega_target, kv_target, ...
'linear' ...
);
F_interp(~valid_mask) = 0;
F_interp(isnan(F_interp)) = 0;
% % t = ((prb_x0(N_elements) + grid_width/2).*sin(angles(iw)))/c0;
% A_w = exp(1i.*omega_target.* t);
% t = reshape(transmit_delay(:,iw),[N_z_axis N_x_axis]);
% A_w = exp(1i.*omega_target.* t');
% 计算权重过滤无效值,大小为(nkx,nkz)
weight = sqrt(k_target.^2 - kv_target.^2);
% invalid_mask = isnan(weight) | ~isreal(weight);
weight(~valid_mask) = 0;
weight(~isreal(weight)) = 0;
G_kz_kx = G_kz_kx + weight.* F_interp;
% G(kx,kz|ku)
G_kz_kx(isnan(G_kz_kx)) = 0;
end
%-----------------------
G_kz_kx = fftshift(G_kz_kx);
image_recon = ifft2(G_kz_kx); % F(kx,kz)→f(x,z)
image_recon = abs(image_recon);
image_recon = image_recon / max(image_recon(:)); % 归一化
image_recon_db = 20*log10(image_recon + eps); % 转换为dB,避免log(0)
figure;
imagesc(1:nkx * dx * 1e3, 1:nkz * dx * 1e3,image_recon);
axis image;
% dbscale = -30;
% clim([dbscale, 0]);
colorbar;
colormap(jet)
grid on;
set(gca, 'GridColor', [0.7, 0.7, 0.7], 'GridLineStyle', ':', 'GridAlpha', 0.5);
% xlabel('Horizontal Position [mm]','fontsize',14)
% ylabel('Depth [mm]','fontsize',14)
title('FW-PWI Image','fontsize',14);
xlabel('x [mm]');
ylabel('z [mm]');