Ультразвуковая визуализация

Question

Здравствуйте, я пытаюсь сделать более быстрый код формирования луча (используя динамическую фокусировку) в Matlab, векторизовав все переменные, которые необходимо обработать. Оригинальный код ниже:

function out = BF(in, maxAprSz, Fs, Fnum)
% The following script take the pre-beamform data and creates beamformed
% image by applying parallel beamformers to the same data. 
%
% input:
%   in:             the first dimension should be scan lines and the second
%                   dimension should be the channel samples
%   maxAprSz:       maximum size of the aperture for beamforming
%   Fs              Acquisition sampling frequency 80*1e6, 40MHz, ... 10*1e6
%   Fnum            desired F number=focal depth/(2*aperture)
% output:
%   out: beamformed rf data

nl = size(in,2);   % number of lines
ns = size(in,1);   % number of samples

pitch = 0.3 * 1e-3;         % spacing between lines/channels in the lateral direction [m]
c = 1540;                   % speed of sound 1540 [m/s]
sampleSpacing = c/Fs/2;     % spacing between samples in the axial direction in milimeter for Fs

out = zeros(nl,ns);

for i = 1:nl    % for each line/channel

  for j=1:ns  % find the value for each sample

    % find the sample location
    Z =  (j ) * sampleSpacing ; %focal depth
    %place of i element in x axis

    % calculate the aperture based on the F number
    fulAprSz = Z/(2*Fnum); 
    hlfAprSz = floor(fulAprSz / pitch);
    % check to make sure we do not go beyond max apr
    if (hlfAprSz > maxAprSz/2)
        hlfAprSz = floor(maxAprSz / 2);
    end

    x = -hlfAprSz: hlfAprSz;    % aperture indices
    fulAprSz = 2*hlfAprSz + 1;

    X = i * pitch;
    X1 = (i + x) * pitch;
    % calc delays based on sample depth and receive aperture

    delays = ( Z + sqrt( Z^2 + (X-X1).^2 ) ) / c * Fs; 
    delays = round(delays); % no interpolation used in this version

    win = hanning(fulAprSz)';   % windowing used for apodization
    chnls = zeros(1,fulAprSz);  % will be filled with proper signal from each channel

    cntr = i;       % center of apreture
    apr = cntr + x; % full aperture indices

    % find the corresponding value from each channel
    for k = 1:fulAprSz  

        chlIndx = apr(k);
        if chlIndx<1, continue, end; 
        if chlIndx>nl, continue, end;

        chlSmpl = delays(k);

        if chlSmpl<1, continue, end;
        if chlSmpl>ns, continue, end;

        chnls(k) = in(chlIndx, chlSmpl); 

    end;

    % apodization : ideally has to be a function of depth
    chnls = win .* chnls;

    % beamforming

    out(i,j) = sum( chnls );
    end;
end

Я думал сделать реализацию с массивами ячеек (потому что fulAprSz изменяется в каждом Z из-за использования Fnum). Я сделал приведенную ниже реализацию, и я рассчитал задержки для одного канала (столбца) для каждого Z (глубина). Эти задержки одинаковы для каждого канала. Я застрял на том, как найти для каждого канала выборки, связанные с этими задержками, а затем суммировать их, чтобы я мог взять выходную диаграмму, сформированную на выборке. В приведенном ниже коде nl=128 (каналы) и ns=2080 (выборки) в качестве примера измерений входных данных. У кого-нибудь есть идеи?

nl=128; %channels
ns=2080; %samples

fs=40e6; %sampling frequency
Fnum=2; %f number

pitch=0.3*1e-3; %space between each channel [m]
c = 1540; %sound of speed [m/s]                 
sampleSpacing = c/fs/2; 

out=zeros(nl,ns);

z=(1:ns)*sampleSpacing; %depth

fulAprSz = z/(2*Fnum); %full aperture size according to fnum
hlfAprSz = floor(fulAprSz / pitch);
fulAprSz = 2*hlfAprSz + 1;

for i=1:ns
    sizeM=fulAprSz(i); %size of each matrix
    win{i}=hanning(sizeM); %apodization matrix for each z
    X1{i}=(-hlfAprSz(i):hlfAprSz(i)); %indices of aperture
    x{i}=X1{i}*pitch; %aperture in meters
    delays{i} = round (( z(i) + sqrt( z(i).^2 + x{i}.^2 ) ) / c * fs); %delays for a single line
end