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