Dear sir
I learnt the 3D FFT Reconstruction For A Planar Sensor Example and want to use the method to reconstructe my real experimental data.However I can not get right result.My code is as follows：
clc
clear all
close all
% start timer
tic;
% nargin=42;
% define defaults
num_req_inputs = 5;
data_order = 'tyz';
interp_method = '*nearest';
plot_recon = 1;
positivity_cond = 1;
c = 1500; % The speed of sound
fs =25; %MHz Sampling rate
dt=1/fs/1000000;
dy=0.00135;
dz=0.0003;
% read data,data order is tyz
load('pa_data_all')
p=pa_data_all(:,:,:);
if strcmp(data_order, 'yzt')
p = permute(p, [3 1 2]);
end
% mirror the time domain data about t = 0 to allow the cosine transform in
% the t direction to be computed using an FFT
p = [flipdim(p, 1); p(2:end, :, :)];
% extract the size of mirrored input data
[Nt, Ny, Nz] = size(p);
% update command line status
disp('Running k-Wave planar reconstruction...');
disp([' grid size: ' num2str((Nt+1)/2) ' by ' num2str(Ny) ' by ' num2str(Nz) ' grid points']);
disp([' interpolation mode: ' interp_method]);
% create a computational grid that is evenly spaced in w, ky, and kz, where
% Nx = Nt and dx = dt*c
kgrid = makeGrid(Nt, dt*c, Ny, dy, Nz, dz);
% from the grid for kx, create a computational grid for w using the
% relation dx = dt*c; this represents the initial sampling of p(w, ky, kz)
w = c*kgrid.kx;
% remap the computational grid for kx onto w using the dispersion
% relation w/c = (kx^2 + ky^2 + kz^2)^1/2. This gives an w grid that is
% evenly spaced in kx. This is used for the interpolation from p(w, ky, kz)
% to p(kx, ky, kz). Only real w is taken to force kx (and thus x) to be
% symmetrical about 0 after the interpolation.
w_new = (c*kgrid.k);
% calculate the scaling factor using the value of kx, where
% kx = sqrt( (w/c).^2 - kgrid.ky.^2 - kgrid.kz.^2 ) and then manually
% replacing the DC value with its limit (otherwise NaN results)
sf = c^2*sqrt( (w/c).^2 - kgrid.ky.^2 - kgrid.kz.^2)./(2*w);
sf(w == 0 & kgrid.ky == 0 & kgrid.kz == 0) = c/2;
% compute the FFT of the input data p(t, y, z) to yield p(w, ky, kz) and
% scale
p = sf.*fftshift(fftn(fftshift(p)));
% remove unused variables
% clear sf;
% exclude the inhomogeneous part of the wave
p(abs(w) < (c*sqrt(kgrid.ky.^2 + kgrid.kz.^2))) = 0;
% compute the interpolation from p(w, ky, kz) to p(kx, ky, kz); for a
% matrix indexed as [M, N, P], the axis variables must be given in the
% order N, M, P
p = interp3(kgrid.ky, w, kgrid.kz, p, kgrid.ky, w_new, kgrid.kz, interp_method);
% remove unused variables
% clear kgrid w;
% set values outside the interpolation range to zero
p(isnan(p)) = 0;
% compute the inverse FFT of p(kx, ky, kz) to yield p(x, y, z)
p = real(ifftshift(ifftn(ifftshift(p))));
% remove the left part of the mirrored data which corresponds to the
% negative part of the mirrored time data
p = p( (Nt + 1)/2:Nt, :, :);
% correct the scaling - the forward FFT is computed with a spacing of dt
% and the reverse requires a spacing of dz = dt*c, the reconstruction
% assumes that p0 is symmetrical about z, and only half the plane collects
% data (first approximation to correcting the limited view problem) (p_zxy)
p = 2*2*p./c;
% enfore positivity condition (p_zxy)
if positivity_cond
disp(' applying positivity condition...');
p(p < 0) = 0;
end
% update command line status
disp([' computation completed in ' scaleTime(toc)]);
% plot the reconstruction
if plot_recon
% allocate axis dimensions
x_axis = [0 (Nt/2)*dt*c];
y_axis = [0 Ny*dy];
z_axis = [0 Nz*dz];
% select suitable axis scaling factor
[x_sc, scale, prefix] = scaleSI(max([(Nt/2)*dt*c, Ny*dy, Nz*dz]));
% select suitable plot scaling factor
plot_scale = max(p(:));
% create the figures
figure;
subplot(2, 2, 1), imagesc(y_axis*scale, x_axis*scale, squeeze(p(:, :, round(end/2))), [-plot_scale, plot_scale]);
axis image;
title('x-y plane');
subplot(2, 2, 2), imagesc(z_axis*scale, x_axis*scale, squeeze(p(:, round(end/2), :)), [-plot_scale, plot_scale]);
axis image;
xlabel(['(All axes in ' prefix 'm)']);
title('x-z plane');
subplot(2, 2, 3), imagesc(z_axis*scale, y_axis*scale, squeeze(p(round(end/2), :, :)), [-plot_scale, plot_scale]);
axis image;
title('y-z plane');
colormap(getColorMap);
end
Thank you Sir.
Sincerely