Введение

Этот пример вводит концепцию фокусируемого луча и показывает, как использовать phased.FocusedSteeringVector сгенерировать необходимые веса элемента для фазированной решетки. Это также показывает, как использовать phased.SphericalWavefrontArrayResponse вычислить ответ массивов массива под данным углом и областью значений. Сначала вы посмотрите на ответ, когда один луч будет сформирован, и исследуйте характеристики фокальной области, затем эмулируйте стратегию набора, обычно используемую в обработке изображений ультразвука, чтобы видеть, как фокусируемый луч появляется в изображении.

Фокусируемый Beamforming

Сферическая задержка-и-сумма фронта импульса

Beamforming под типичным предположением далекого поля позволяет моделировать фронт импульса сигнала как плоскость, которая распространяет вдоль ее нормального направления. Когда сигнал является инцидентом на массиве, фронт импульса пересекает отдельные элементы с относительной задержкой, которая пропорциональна расстоянию элемента вдоль направления распространения волны. В соответствии с этой моделью, задержка (или сдвиг фазы, в узкополосном случае) может быть применена к каждому элементу, таким образом, что выходные параметры через элементы (получают ли на передаче или) имеют конструктивную фазу, когда когерентно суммировано.

Без предположения далекого поля фронт импульса, происходящий от точечного источника, моделируется как сферическая поверхность, сосредоточенная в том Соре. Этот фронт импульса пересекает элементы фазированной решетки с относительной задержкой, продиктованной гиперболической областью значений через колинейные элементы. В то время как эта модель просто добавила бы ненужную стоимость расчета для генерации далекой диаграммы направленности по напряжённости поля, сферическая модель фронта импульса необходима, чтобы изучить почти поле beamforming и фокусировку.

Функциональный helperPlotULAWavefronts обеспечивается, чтобы продемонстрировать различие в задержках элемента и форме фронта импульса между управляемым лучом с и без фокусировки. Задержки элемента относительно центра массивов.

figure
helperPlotULAWavefronts(14,1e6,1540,20,inf)
title('Steered Wavefront')

figure
helperPlotULAWavefronts(14,1e6,1540,20,0.01)
title('Steered and Focused Wavefront')

Фокальная область и Близкий/Далекий Контур

В далеком поле управляемый массив имеет четко определенный шаблон на угловом пробеле, который не зависит от области значений. В близком поле управляемый (но нефокусируемый) массив не имеет никакой discernabe структуры лепестка вообще. Благодаря нелинейному фазовому соотношению между элементами ответ, который имеет равную величину во всех областях значений в данном направлении, не возможен. Вместо этого путем фокусировки луча можно получить небольшую область, ограниченную и в углу и в области значений, в которой ответ напоминает ответ далекого поля, известный как фокальную область. Угловым положением фокальной области можно управлять так же легко как с лучом далекого поля. Положение фокальной области в области значений определяется фокальной областью значений и глубиной резкости (DoF), которая является степенью области значений области.

Рисунки ниже показывают прогрессию формы луча с областью значений для фокусируемых универсальных и прямоугольных массивов. Функциональный helperPlotResponseSlices обеспечивается, чтобы продемонстрировать, как сгенерировать этот тип фигуры.

В фокальной области значений наш луч на угловом пробеле тесно напоминает луч далекой диаграммы направленности по напряжённости поля.

Близкий/далекий контур, подобный в концепции к другим near-field/far-field контурам, является контуром мимо, какая фокусировка не возможна. С другой стороны управляемый, но нефокусируемый луч не возможен на близкой стороне контура. Рисунок ниже демонстрирует этот факт. Заметьте, что управляемый луч только начинает формироваться на противоположной стороне, и фокусируемый луч только фокусируется на близкой стороне.

Область значений, где этот контур может быть найден, хорошо зарегистрирована в медицинскую литературу обработки изображений как ${\mathit{L}}_{\mathrm{array}}^2/4\lambda$, где ${\mathit{L}}_{\mathrm{array}}$ длина массива. Для простой универсальной линейной матрицы с $\mathit{N}$ критически распределенные элементы, это может быть описано как ${\mathit{N}}^2\lambda /16$.

Генерация ответа для одного луча

Это разделяет, показывает, как сгенерировать и построить один луч. Запустите путем подготовки фазированной решетки и других Системных объектов. Общие медицинские системы ультразвука обработки изображений действуют между 2 и 20 МГц, и обычно принятая средняя скорость распространения звука в мягкой ткани составляет 1 540 м/с.

rng('default')
freq = 4e6;
c = 1540;
lambda = c/freq;

Преобразователи ультразвука прибывают в большое разнообразие топологии, которой удовлетворяют для определенных действий. Этот пример просто использует универсальную линейную матрицу с критическим интервалом элемента.

numElems = 256;
elemSpacing = lambda/2;
array = phased.ULA(numElems,elemSpacing);

Системные объекты phased.FocusedSteeringVector и phased.SphericalWavefrontArrayResponse используются в тандеме, чтобы сгенерировать управляемые и фокусируемые веса элемента и вычислить ответ по некоторой области. Передайте массив и спецификацию скорости распространения сверху конструкторам. Для ответа массивов также включите входной порт весов.

SV = phased.FocusedSteeringVector('SensorArray',array,'PropagationSpeed',c);
AR = phased.SphericalWavefrontArrayResponse('SensorArray',array,'PropagationSpeed',c,'WeightsInputPort',true);

Теперь у вас есть минимальная настройка, требуемая сформировать и смотреть фокусируемый луч. Использование 10 степеней держится в азимуте и фокальной области значений 40 мм.

azSteer = 10;
focalRange = 0.04;

Используйте область, которая запускается перед массивом и расширяет к дважды фокальной области значений и покрывает длину массива в боковом направлении. Элементы массива простираются вдоль оси Y и массива, нормальное направление является +X.

arrayLength = numElems*elemSpacing;
x = linspace(1e-3,2*focalRange,200);
y = linspace(-arrayLength/2,arrayLength/2,200);

Преобразуйте в сферические координаты для входа к держащемуся вектору и расчету ответа.

[az,el,rng] = cart2sph(x,y',0);
ang = rad2deg([az(:) el(:)]');
rng = rng(:)';

Теперь вычислите ответ фокусируемого массива по заданной области.

weights = SV(freq,[azSteer;0],focalRange);
beam = AR(freq,ang,rng,weights);

Измените, нормируйте и используйте логарифмическую шкалу.

beam = reshape(beam,numel(y),numel(x));
beam = beam/max(abs(beam(:)));
beam = mag2db(abs(beam));

Используйте обеспеченный функциональный helperPlotResponse построить результат. Положения элементов массива обозначаются красными маркерами.

figure
helperPlotResponse(beam,x,y,array)
title('Single Steered and Focused Beam')

Фокальная область ясно отображается в заданной области и углу. Как с ответом далекого поля (шаблон), величина сферического ответа фронта импульса учитывает различную фазу через элементы, но, в отличие от ответа далекого поля, сферический ответ фронта импульса также включает эффект различной потери распространения свободного пространства через элементы. Это - по существу небольшая амплитудная модуляция через элементы, эффект которых отображается в боковых лепестках луча. Плоское усиление, равное фокальной области значений, применяется к каждому элементу, так, чтобы полное амплитудное взвешивание на элементе было отношением расстояния между точкой ответа, и центр массивов к расстоянию от ответа указывают на отдельный элемент: ${\mathit{R}}_{\mathrm{resp}}/{\mathit{R}}_{\mathrm{elem}}$. Например, вклад в луч суммы от элемента, расположенного в начале координат, имел бы модульную величину.

Фокальная область

Для фокальной области, которая будет ограничена в области значений, фокальная область значений должна быть малой достаточно, что весь регион (со степенью, определенной степенью свободы), должен быть ближе, чем близкий/далекий контур. Степень свободы может быть описана в зависимости от связанного количества, обычно замечаемого в оптике, известной как F-номер, который является отношением фокальной области значений к длине массива: $\mathit{F}={\mathit{R}}_{\mathrm{foc}}/{\mathit{L}}_{\mathrm{array}}$. От этого количества хорошая оценка степени свободы ${\mathit{d}}_{\mathit{F}}=\mathrm{7}.1\lambda {\mathit{F}}^2$. Для фиксированного массива и частоты, степень свободы увеличивается с квадратом фокальной области значений.

Этот раздел смотрит на то, как степень свободы изменяется с фокальной областью значений. Сгенерируйте ответ на интервале фокальных областей значений и исследуйте изменяющуюся степень свободы. Функциональный helperPlotBeamMarkers отобразит индикацию относительно степени свободы для каждой фокальной области значений. Фокальная область является примерно эллиптической, и не сосредоточена на фокальной области значений. Другое эмпирическое правило может использоваться, чтобы найти смещение от фокальной области значений до центра области: ${\mathit{d}}_{\mathit{F}}/4$.

nearField = arrayLength^2/(4*lambda); % Near/far boundary range

x = linspace(1e-3,nearField*1.2,100);
y = linspace(-0.01,0.01,100);

coeff = 0.1:0.1:0.9; % Proportion of near/far boundary
focalRanges = coeff*nearField;

fnum = focalRanges/arrayLength;
dof = 7.1*lambda*fnum.^2;
focalRegionCenter = focalRanges + dof/4;

figure
for ind = 1:numel(focalRanges)
    beam = helperMakeSingleBeam( SV,AR,freq,0,focalRanges(ind),x,y );
    helperPlotResponse(beam,x,y)
    caxis([-6 0]) % only view the greatest 6 dB of the beam
    helperPlotBeamMarkers(focalRanges(ind),focalRegionCenter(ind),nearField,dof(ind),0.001)
    title(sprintf('Focal Range = %d%% of Near/Far Boundary',round(coeff(ind)*100)))
    drawnow
end

К тому времени, когда фокальная область значений увеличилась до 60% близкого/далекого контура, эксцентриситет фокальной области стал довольно большим. На 80% фокальная область пересекла близкий/далекий контур и по существу стала нефокусируемым лучом.

Ширина луча в фокальной области, как с ответом далекого поля, может все еще быть примерно вычислена с обычным $\lambda /{\mathit{L}}_{\mathrm{array}}$, отношение длины волны к длине массива. Таким образом ширина (в боковом направлении) фокальной области может быть аппроксимирована ${\mathit{R}}_{\mathrm{foc}}\lambda /{\mathit{L}}_{\mathrm{array}}$.

Однострочный захват с линейной переменой подрешетки

A-сканы и формирование изображений

В отличие от радиолокационных систем, где направление прибытия может быть оценено эффективно через beamforming сигналов далекого поля, почти поле beamforming и требования задержки систем ультразвука часто требуют более простой стратегии определить местоположение источника возвращенной энергии. Общий класс стратегии набора включает несколько A-сканов, профилей отражающей способности вдоль данного линейного сегмента. Размещение этих линий в сформированном изображении определяется просто известным положением луча.

Сформировать прямоугольное изображение, как сделан в этом примере, lineary перемена подрешетки, может использоваться. С этим методом подмножество элементов массива (подрешетка) используется для каждого импульса, без регулирования, чтобы получить профиль области значений, происходящий в центре подрешетки и расширяющий в осевом направлении. Последовательные линии формируются путем перемещения выбора подрешетки к различному набору элементов.

Выбор фокальной области значений может быть рассмотрен отдельно от обращения луча. Некоторые системы могут сохранить фиксированную фокальную область значений или использовать динамическую фокусировку на, получают наряду с аподизацией или работой с окнами, чтобы сгенерировать нечеткое изображение, которое покрывает большую часть почти полевой области. В этом примере фокальная область значений сохранена фиксированной, в то время как боковое положение подрешетки варьируется.

Класс помощника helperSubarray обеспечивается, чтобы эмулировать выбор подрешетки и упростить цикл симуляции. Этот класс отслеживает, которых элементы принадлежат текущей подрешетке, и обрабатывает необходимые преобразования между системами координат подрешетки и глобальной переменной.

Отобразите формирование

Для того, чтобы продемонстрировать эффекты фокусируемого beamforming, этот пример использует простую стратегию формирования изображений, которая создает каждый профиль области значений путем квантования и накопления ответа в каждом импульсе, затем вставки того профиля в завершенное изображение на основе местоположения подрешетки. Это симулирует идеальный импульс дельты в свободном пространстве, которое допускает сравнение бокового разрешения внутри и снаружи фокальной области. Этот метод игнорирует эффекты выбора формы волны и многопутевых отражений, и обрабатывает рассеиватели как совершенную точку изотропные отражатели.

Симуляция

Те же системные параметры будут использоваться в качестве в предыдущем разделе. Каждая подрешетка будет состоять из 64 элементов. Подрешетка запускается в конце массива и смещена одним элементом на каждом импульсе.

numSubElems = 64;
subarray = helperSubarray(array,numSubElems);

Если подрешетка будет смещена одним элементом на каждом импульсе, и все непрерывные подрешетки покрыты, общее количество импульсов будет

numPulses = numElems - numSubElems + 1

numPulses = 193

После анализа в предыдущем разделе проверяйте фокальные параметры области на эту подрешетку. Получите длину подрешетки, близкий/далекий контур, степень свободы и боковую ширину фокальной области.

subarrayLength = numSubElems*elemSpacing;
nearFieldSub = subarrayLength^2/(4*lambda); % Near/far boundary for the subarray
fnumSub = focalRange/subarrayLength;
dofSub = 7.1*lambda*fnumSub.^2

dofSub = 0.0288

widthSub = lambda/subarrayLength*focalRange

widthSub = 0.0013

Луч подрешетки имеет степень свободы приблизительно 28,8 мм, и боковая ширина фокальной области составляет приблизительно 1,3 мм. Проверяйте, что фокальная область хорошо в близком/далеком контуре

boundedFocalRegion = focalRange + dofSub < nearFieldSub

boundedFocalRegion = logical
   1

Чтобы продемонстрировать первичную полноценность фокусируемого луча, уменьшенной ширины луча в фокальной области, используют несколько параллельных линий рассеивателей вдоль осевого направления (линии глубины). Задайте желаемую макс. глубину рассеивателей и используйте обеспеченную функцию помощника, helperGetResponsePoints, получить положения рассеивателя. Позвольте всем рассеивателям иметь отражающую способность с модульной амплитудой. Положения рассеивателя встревожены, чтобы избежать artifacting из-за симметрии.

Для бокового интервала линий используйте 4 раза расчетную фокальную ширину области. Поскольку подрешетка просто смещена одним элементом за один раз, который является более коротким расстоянием, чем наша ширина луча в фокальной области, возвратитесь из каждой линии рассеивателей, обнаруживается больше чем в одной строке изображения, заставляя их, казаться, иметь большую ширину.

maxDepth = nearFieldSub;
lineSpacing = 4*widthSub;
[sx,sy] = helperGetResponsePoints(maxDepth,arrayLength,lambda,lineSpacing);

Чтобы визуализировать сцену, постройте положения рассеивателя наряду с элементами массива.

figure
plot(subarray)
hold on
plot(sx,sy,'.')
hold off
legend('Array Elements','Initial Subarray','Scatterers')
xlabel('Axial Distance')
ylabel('Lateral Position')

Чтобы сформировать профиль области значений и получить эффекты интерференции от бокового лепестка возвращаются, параметры выборки области значений должны быть заданы. Современные системы ультразвука используют относительно высокую частоту дискретизации того же порядка как центральная частота переданной формы волны. Используйте размер интервала области значений 1 мм, соответствуя частоте дискретизации приблизительно 1,5 МГц. Предоставленный помощник функционирует helperFormRangeProfile используется, чтобы сформировать профиль диапазона из данных об ответе массивов.

rangeBinSize = 1e-3;
Fs = c/rangeBinSize;
rangeBins = 0:rangeBinSize:maxDepth;
numRangeSamples = numel(rangeBins);

Получите положения рассеивателя в сферических координатах (углы в градусах) для входа к SphericalWavefrontArrayResponse.

[respAng,~,respRng] = cart2sph(sx,sy,0);
respAng = rad2deg(respAng(:)');
respRng = respRng(:)';

Каждый цикл симуляции включает вычислительные веса для текущей подрешетки, вычисляя ответ, формируя профиль области значений, и формируя изображение линию за линией. Мы также применим зависимое областью значений усиление для универсальной яркости (см. helperFormRangeProfile).

im = zeros(numPulses,numRangeSamples);
center = zeros(3,numPulses);

for pulse = 1:numPulses
    
    center(:,pulse) = subarray.center; % Center position of our current subarray
    [focAzGlobal,~,focRngGlobal] = subarray.localToGlobalSph( 0,0,focalRange ); % Angle and range to current focal point
    
    weights = SV(freq,[focAzGlobal;0],focRngGlobal);
    weights(~subarray.selection) = 0; % Zero-out weights for unused elements

    resp = AR(freq,respAng,respRng,weights);
    resp = resp./respRng(:); % Undo normalization to use actual propagation loss

    im(pulse,:) = helperFormRangeProfile(resp,sx,sy,center(:,pulse),rangeBins); % Add line to image

    if pulse < numPulses
        subarray.shift(1) % Shift subarray if there are pulses remaining
    end

end

Нормируйте и постройте изображение, наряду с наложением положений рассеивателя.

im = im/max(abs(im(:)));
figure
subplot(1,2,1)
helperPlotResponse(mag2db(abs(im)),rangeBins,center(2,:))
title('Linear Subarray Shift Image')
subplot(1,2,2)
helperPlotResponse(mag2db(abs(im)),rangeBins,center(2,:))
hold on
plot(sx,sy,'.r','markersize',1)
hold off
title('Scatterer Overlay')
set(gcf,'Position',get(gcf,'Position')+[0 0 560 0]);

Линии глубины разрешимы о фокальной области значений, на интервале области значений примерно равняются вычисленной степени свободы для луча подрешетки. Далеко от центра, параллельные линии рассеивателей быстро становятся неразличимыми благодаря широкому угловому распространению луча за пределами фокальной области.

Обратите внимание на то, что, хотя все рассеиватели были использованной для расчета энергией, возвращаются, не, все линии отображаются из-за усечения между полным размером массивов и степенью позиций центра подрешетки. Более широкое общее поле зрения было бы доступно с более короткой подрешеткой, за счет уменьшаемого бокового разрешения.

Заключение

Этот пример ввел два Системных объекта для вычисления фокусируемых весов и для вычисления ответа "не далекое поле" массива со сферическими фронтами импульса. Это показало, как исследовать некоторые основные характеристики фокусируемого луча, и как сгенерировать основное изображение, чтобы визуализировать эффект фокусируемого луча на боковом разрешении с линейным методом сбора сдвига подрешетки.

Ссылки

[1] Деми, Libertario. “Практическое Руководство по Формированию Луча Ультразвука: Диаграмма направленности и Анализ Реконструкции Изображений”. Прикладные науки 8, № 9 (3 сентября 2018): 1544. https://doi.org/10.3390/app8091544

[2] Ramm, О. Т. Фон и С. В. Смит. “Излучите Регулирование с Линейными матрицами”. Транзакции IEEE на Биоинженерии BME-30, № 8 (август 1983): 438–52.

Функции помощника

helperMakeSingleBeam

function [beam,x,y] = helperMakeSingleBeam( SV,AR,freq,azSteer,focalRange,x,y )
% Get the element weighting for a single beam and compute the response

weights = SV(freq,[azSteer;0],focalRange);

% Domain of response
if nargin < 6
    x = linspace(1e-3,2*focalRange,200);
end
if nargin < 7
    pos = SV.SensorArray.getElementPosition;
    halfy = max(pos(2,:))*1.2;
    y = linspace(-halfy,halfy,200);
end
[az,el,rng] = cart2sph(x,y',0);
ang = rad2deg([az(:) el(:)]');
rng = rng(:)';

% Generate response
beam = AR(freq,ang,rng,weights);
beam = reshape(beam,numel(y),numel(x));

% Normalize and use log scale
beam = beam/max(abs(beam(:)));
beam = mag2db(abs(beam));

end

helperPlotResponse

function helperPlotResponse(R,x,y,array)
% Plot the response, R, on the domain defined by x and y.

imagesc(x,y,R)
set(gca,'ydir','normal')
caxis([-32 0])
xlabel('Axial Distance')
ylabel('Lateral Position')

if nargin > 3
    pos = getElementPosition(array);
    hold on
    h = plot(pos(1,:),pos(2,:),'.r');
    hold off
    xl = xlim;
    xlim([min(pos(1,:)) xl(2)])
    legend(h,'Array Elements','Location','southeast','AutoUpdate','off')
end

end

helperPlotBeamMarkers

function helperPlotBeamMarkers(focalRange,center,nearField,dof,offset)
% Put informative markers on a beam plot

line([center-dof/2 center+dof/2],[offset offset],'color','white')
line([center-dof/2 center-dof/2],[0 offset],'color','white')
line([center+dof/2 center+dof/2],[0 offset],'color','white')

line([focalRange focalRange],ylim,'color','red')
line([nearField nearField],ylim,'color','cyan')

end

helperGetResponsePoints

function [sx,sy] = helperGetResponsePoints( maxDepth,arrayLength,lambda,dy )
% Make parallel lines of scatterers along X

sx = linspace(0.001,maxDepth,400);
sy = -arrayLength/2:dy:arrayLength/2;

[sx,sy] = meshgrid(sx,sy);
sx = sx(:);
sy = sy(:);

sx = sx + (rand(size(sx))-1/2)*lambda;
sy = sy + (rand(size(sy))-1/2)*lambda;

end

helperFormRangeProfile

function rangeProf = helperFormRangeProfile(resp,sx,sy,center,rangeBins)
% This helper function quantizes a response in range, coherently
% accumulates the return in each bin, and applies amplitude weighting per
% range bin

rangeBinSize = rangeBins(2) - rangeBins(1);
numRangeSamples = numel(rangeBins);

% Range of scatterers relative to subarray center
scatRngRel = sqrt((center(1)-sx).^2 + (center(2)-sy).^2);

% Quantize scatterer ranges into fast-time sampling vector
scatRidx = 1 + floor(scatRngRel/rangeBinSize);

% Only keep samples below max depth
I = scatRidx <= numRangeSamples;
scatRidx = scatRidx(I);
resp = resp(I);

% Accumulate return into fast-time sampling grid
rangeProf = accumarray(scatRidx,resp,[numRangeSamples 1]);
rangeProf = rangeProf';

% Apply range-dependent gain
rangeProf = rangeProf.*rangeBins;

end

helperPlotULAWavefronts

function helperPlotULAWavefronts( numElems,f,c,az,r )
% Plot the wavefronts for the given ULA with ArrayAxis 'y', for the given
% azimuth angle and focal range.
%
% For the far-field wavefront, use inf for focal range

lambda = c/f;
array = phased.ULA(numElems,lambda/2);
pos = getElementPosition(array);
arrayLength = max(pos(2,:)) - min(pos(2,:));

% get relative path lengths
if isinf(r)
    L = phased.internal.elemdelay(pos,c,[az;0])*c;
else
    L = phased.internal.sphericalelemdelay(pos,c,[az;0],r)*c;
end

% plot element positions
plot(pos(1,:),pos(2,:),'oblue');
hold on;

% if near field, plot source
if ~isinf(r)
    [src(1,1),src(2,1),src(3,1)] = sph2cart(az*pi/180,0,r);
    plot(src(1),src(2),'*r','markersize',10);
end

% wavefront marker width
s = lambda/6;

% far field prop path
if isinf(r)
    [los(1,1),los(2,1),los(3,1)] = sph2cart(az*pi/180,0,1);
end

for ind = 1:size(pos,2) % for each element
   
    p = pos(:,ind);

    if isinf(r)
        src = p + los*arrayLength;
    end
    
    path = src - p;
    path = path/norm(path);
    
    wp = p + path*L(ind); % position of wavefront
    
    % prop path
    line([wp(1) src(1)],[wp(2) src(2)],'color','black','linestyle','--');
    
    % wavefront marker
    u = cross([0;0;1],path)*s;
    line([wp(1)-u(1) wp(1)+u(1)],[wp(2)-u(2) wp(2)+u(2)],'color','magenta');
    
end

hold off;
grid on;
axis equal;

if isinf(r)
    legend('Elements','Prop Path','Wavefronts','Location','SouthEast');
else    
    legend('Elements','Focal Point','Prop Path','Wavefronts','Location','SouthEast');
end

end

helperPlotResponseSlices

function helperPlotResponseSlices
% Demonstrates how to visualize range slices of a spherical wavefront response

f = 2e6;
c = 1540;
lambda = freq2wavelen(f,c);

array = phased.URA([32 32],lambda/2);
elemPos = array.getElementPosition;

focalRange = 0.03;
sampleRanges = .01:.01:.05;

% domain of each slice
azSteer = -20;
elSteer = 20;
az = azSteer + (-30:.1:30);
el = elSteer + (-30:.1:30);
[az,el] = meshgrid(az,el);
ang = [az(:) el(:)]';
[x,y,z] = sph2cart(az*pi/180,el*pi/180,1);

SV = phased.FocusedSteeringVector('SensorArray',array,'PropagationSpeed',c);
AR = phased.SphericalWavefrontArrayResponse('SensorArray',array,'PropagationSpeed',c,'WeightsInputPort',true);

w = SV(f,[azSteer;elSteer],focalRange);

for ind = 1:numel(sampleRanges)
    resp = AR(f,ang,sampleRanges(ind),w);
    resp = resp / array.getNumElements;
    resp = reshape(resp,size(x));
    alpha = 1 - (abs(sampleRanges(ind) - focalRange)/(max(sampleRanges)-min(sampleRanges))); % transparency

    surf(x*sampleRanges(ind),y*sampleRanges(ind),z*sampleRanges(ind),mag2db(abs(resp)),'FaceAlpha',alpha)
    hold on
    shading flat
end

% plot element positions and boresight vector
caxis([-32 0])
plot3(elemPos(1,:),elemPos(2,:),elemPos(3,:),'.black','markersize',4);
[b(1),b(2),b(3)] = sph2cart(azSteer*pi/180,elSteer*pi/180,max(sampleRanges)*1.2);
quiver3(0,0,0,b(1),b(2),b(3),'black','autoscale','off')
axis equal
axis off
hold off
set(gca,'view',[-70 22])

end