Model Hybrid Beamforming with CDL Channel

Since R2025a

This example shows how to apply digital and analog beamforming to a 5G New Radio (NR) physical downlink shared channel (PDSCH) link by using a CDL channel.

Introduction

Fully digital beamforming techniques equip one radio frequency (RF) chain for each antenna element. This architecture presents challenges in terms of power consumption, cost, and hardware complexity, especially at millimeter-wave frequencies. Hybrid beamforming combines analog and digital processing to optimize performance while reducing the number of RF chains, power usage, and hardware requirements. This approach enables flexible beam management and spatial multiplexing, where multiple data streams are transmitted simultaneously. This example illustrates how to combine digital and analog beamforming techniques in a 5G NR PDSCH transmission.

Two widely considered approaches are the fully and partially connected architectures. In the fully connected architecture, each RF chain connects to all elements of an antenna array. In the partially connected architecture, each RF chain connects to a subset of antenna elements. This example demonstrates a partially connected architecture by using fixed rectangular subarrays of equal size.

Simulation Parameters

Set the signal-to-noise ratio (SNR) for the simulation. The SNR for each layer and resource element (RE) accounts for the signal and noise across all antennas. For an explanation of the SNR definition that this example uses, see the SNR Definition Used in Link Simulations example.

simParameters = struct();
simParameters.SNRIn = 10; % dB

Channel Estimator Configuration

The logical variable PerfectChannelEstimator controls the channel estimation behavior. When you set this variable to true, the simulation uses perfect channel estimation. When you set it to false, the simulation uses practical channel estimation based on the values of the received PDSCH demodulation reference signal (DM-RS).

simParameters.PerfectChannelEstimator = true;

Carrier and PDSCH Configuration

Set a 10 MHz bandwidth carrier with a 15 kHz subcarrier spacing.

% SCS carrier parameters
simParameters.Carrier = nrCarrierConfig;         % Carrier resource grid configuration
simParameters.Carrier.NSizeGrid = 52;            % Bandwidth in number of resource blocks
simParameters.Carrier.SubcarrierSpacing = 15;    % 15, 30, 60, 120 (kHz)

Specify a multi-layer full-band slot-wise PDSCH transmission.

simParameters.PDSCH = nrPDSCHConfig;

% Define PDSCH time allocation in a slot (Mapping Type A)
simParameters.PDSCH.SymbolAllocation = [0 simParameters.Carrier.SymbolsPerSlot]; % Starting symbol and number of symbols of each PDSCH allocation

% Define PDSCH frequency resource allocation per slot to be full grid 
simParameters.PDSCH.PRBSet = 0:simParameters.Carrier.NSizeGrid-1;

% Define the number of transmission layers to be used
simParameters.PDSCH.NumLayers = 2;

Antenna Array Configuration

Configure rectangular antenna arrays by specifying their size and number of polarizations. In this example, the antenna array consists of nonoverlapping rectangular subarrays of equal size. You can configure the number of vertical and horizontal elements of the antenna array (${\mathit{N}}_{\mathit{V}}$, ${\mathit{N}}_{\mathit{H}}$) and of the subarrays ($$N_V^S$$, $$N_H^S$$). The antenna array and subarray sizes must be compatible to fit an integer number of subarrays in the array, that is, $$N_V/N_V^S$$ and $$N_H/N_H^S$$ must be integers. This diagram illustrates the partitioning of a 4-by-4 antenna array into 2-by-2 subarrays:

simParameters.TransmitAntennaArray.Size             = [4 4];     % Number of vertical and horizontal antenna elements in array
simParameters.TransmitAntennaArray.SubarraySize     = [2 2];     % Number of vertical and horizontal elements in subarray
simParameters.TransmitAntennaArray.ElementSpacing   = [0.5 0.5]; % Vertical and horizontal distance between antenna elements (wavelengths)
simParameters.TransmitAntennaArray.NumPolarizations = 2;         % Number of polarizations

simParameters.ReceiveAntennaArray.Size              = [2 1];     % Number of vertical and horizontal antenna elements
simParameters.ReceiveAntennaArray.ElementSpacing    = [0.5 0.5]; % Vertical and horizontal distance between antenna elements (wavelengths)
simParameters.ReceiveAntennaArray.NumPolarizations  = 2;         % Number of polarizations

This diagram illustrates three different ways of partitioning a 4-by-4 antenna array into 4 vertical subarrays of size 4-by-1 (left), 8 vertical subarrays of size 2-by-1 (center), and 2 horizontal subarrays of size 2-by-4 elements (right):

To configure single-input single-output (SISO) transmissions with single-polarized antenna elements, set the antenna array and subarray sizes to [1 1] and the number of polarizations to 1.

RF Connections to Antenna Elements

Configure how each RF chain connects to each antenna element within a subarray. Each antenna subarray has two polarizations, and each polarization is fully connected to a single RF chain. This example calculates the connections between the RF chains and their corresponding antenna elements for this architecture using an array of subarrays. For other architectures, you can specify your own connections by using a column vector or matrix.

rfc = RFChainConnections(simParameters.TransmitAntennaArray);

To specify the RF connections to antenna elements, use the value of the connections vector rfc(:) at the ith row to indicate the RF chain index for the ith antenna element within the antenna array. For elements connected to multiple RF chains, you can specify a matrix where the ith row contains the indices of the RF chains that the ith element is connected to. If you connect multiple RF chains to the same antenna element, the transmitter adds the signals before transmission. This diagram describes the subarray connections of a 4-by-2 antenna array partitioned into 2 subarrays of size 2-by-2 antenna elements:

Display the number of RF chains connected to the transmit antenna array.

simParameters.TransmitAntennaArray.RFChainConnections = rfc;
Nrfc = length(unique(rfc));
disp("Number of RF chains: " + Nrfc)

Number of RF chains: 8

simParameters.TransmitAntennaArray.NumRFChains = Nrfc;

Plot how each RF chain connects to each subarray of the transmit antenna array. Each color represents an RF chain that connects to a subarray of antenna elements. The plot displays polarizations in separate figures.

plotRFConnectionsToSubarrays(rfc,simParameters.TransmitAntennaArray)

Obtain and display the overall number of antennas in each array.

simParameters.NTxAnts = numAntennaElements(simParameters.TransmitAntennaArray);
simParameters.NRxAnts = numAntennaElements(simParameters.ReceiveAntennaArray);
disp("Number of transmit antennas: " + simParameters.NTxAnts)

Number of transmit antennas: 32

disp("Number of receive antennas: " + simParameters.NRxAnts)

Number of receive antennas: 4

Propagation Channel Configuration

Configure the delay profile, delay spread, and maximum Doppler shift of the CDL propagation channel for the simulation. You can configure the antenna array geometry in the Antenna Array Configuration section.

simParameters.DelayProfile = 'CDL-C';   % 'CDL-A',...,'CDL-E'
simParameters.DelaySpread = 300e-9;     % s
simParameters.MaximumDopplerShift = 5;  % Hz

simParameters.Channel = createChannel(simParameters);

validateNumLayers(simParameters)

Digital Precoder and Analog Beamformer

In a 5G NR system, the gNB uses channel state information (CSI) from the user equipment (UE) or uplink channel estimates to select suitable digital precoders and analog beamformers. For UE-reported CSI feedback, the gNB can beamform multiple CSI-RS processes independently with a predefined codebook of analog beamformers. Then, the gNB can use the CSI feedback to select appropriate analog beamformers in subsequent data transmissions. Instead of CSI-based precoding and beamforming, this example uses a DFT-based digital precoder. For more information on CSI-based and reciprocity-based digital precoding, see the NR PDSCH Throughput Using Channel State Information Feedback and TDD Reciprocity-Based PDSCH MU-MIMO Using SRS examples.

% Define digital precoder
Nrfc = simParameters.TransmitAntennaArray.NumRFChains;
W = dftmtx(max(simParameters.PDSCH.NumLayers,Nrfc));
W = W(1:simParameters.PDSCH.NumLayers,1:Nrfc);

Specify the analog beamformer. You can configure the azimuth and downtilt angles in degrees of the main beam for each subarray. The azimuth and downtilt angles are relative to the boresight of the array. If you define scalar azimuth and downtilt angles, you use the same beamformer for all subarrays.

subarrayDims = simParameters.TransmitAntennaArray.SubarraySize;
elementSpacing = simParameters.TransmitAntennaArray.ElementSpacing;
numRFChains = simParameters.TransmitAntennaArray.NumRFChains;

azimuth  = linspace(-60,60,numRFChains); % deg
downTilt = 5;                            % deg

F = analogBeamformer(subarrayDims,elementSpacing,numRFChains,azimuth,downTilt);

Display the radiation pattern associated to the first four beamformers.

plotBeams(subarrayDims,elementSpacing,F,1:min(4,numRFChains));

Transmit and Receive PDSCH with Hybrid Beamforming

Apply the digital precoder and analog beamformer to a PDSCH transmission, pass the CP-OFDM waveform through the CDL channel, and recover the transmitted data.

% Reset the random number generator for repeatability
rng("default");

% Extract carrier and PDSCH configuration parameters
carrier = simParameters.Carrier;
pdsch = simParameters.PDSCH;

% Set up propagation channel and receiver
[channel,maxChDelay] = setupChannel(simParameters);
[N0,noiseEst] = setupReceiver(simParameters,channel);

% Calculate the RE capacity for PDSCH allocation
[pdschIndices,pdschIndicesInfo] = nrPDSCHIndices(carrier,pdsch);

% New block of data
data = randi([0 1],pdschIndicesInfo.G,1);

% Create an OFDM resource grid for a slot
dlGrid = nrResourceGrid(carrier,simParameters.TransmitAntennaArray.NumRFChains);

% PDSCH modulation, digital precoding and mapping
pdschSymbols = nrPDSCH(carrier,pdsch,data);
[pdschAntSymbols,pdschAntIndices] = nrPDSCHPrecode(carrier,pdschSymbols,pdschIndices,W);
dlGrid(pdschAntIndices) = pdschAntSymbols;

% PDSCH DM-RS digital precoding and mapping
dmrsSymbols = nrPDSCHDMRS(carrier,pdsch);
dmrsIndices = nrPDSCHDMRSIndices(carrier,pdsch);
[dmrsAntSymbols,dmrsAntIndices] = nrPDSCHPrecode(carrier,dmrsSymbols,dmrsIndices,W);
dlGrid(dmrsAntIndices) = dmrsAntSymbols;

% OFDM modulation
txWaveform = nrOFDMModulate(carrier,dlGrid);

% Analog beamforming of time-domain transmit waveform
txWaveform = analogBeamformWaveform(txWaveform,F,rfc,simParameters.NTxAnts);

% Pass waveform through propagation channel
txWaveform = [txWaveform; zeros(maxChDelay,size(txWaveform,2))];
[rxWaveform,ofdmResponse,timingOffset] = channel(txWaveform,carrier);

% Add AWGN to the received time-domain waveform
noise = N0*randn(size(rxWaveform),'like',1i);
rxWaveform = rxWaveform + noise;

% Synchronize the received waveform
rxWaveform = rxWaveform(1+timingOffset:end,:);

% OFDM demodulate the received waveform
rxGrid = nrOFDMDemodulate(carrier,rxWaveform);

if simParameters.PerfectChannelEstimator
    % Combine perfect channel estimate to account for analog beamformers
    Hest = combineChannelEstimate(ofdmResponse,F,rfc,simParameters.NTxAnts);

    % Get PDSCH resource elements from the received grid and channel
    % estimate
    [pdschRx,pdschHest,~,pdschHestIndices] = nrExtractResources(pdschIndices,rxGrid,Hest);

    % Apply digital precoding to channel estimate
    pdschHest = nrPDSCHPrecode(carrier,pdschHest,pdschHestIndices,W.');
else
    % Practical channel estimation between the received grid and
    % each transmission layer, using the PDSCH DM-RS for each
    % layer. This channel estimate includes the effect of
    % digital precoding and analog beamforming
    [Hest,noiseEst] = nrChannelEstimate(carrier,rxGrid,dmrsIndices,dmrsSymbols,'CDMLengths',pdsch.DMRS.CDMLengths);

    % Get PDSCH resource elements from the received grid and channel
    % estimate
    [pdschRx,pdschHest] = nrExtractResources(pdschIndices,rxGrid,Hest);
end

% Equalize received PDSCH symbols
[pdschEq,eqCSIScaling] = nrEqualizeMMSE(pdschRx,pdschHest,noiseEst);

% Decode PDSCH physical channel
[llr,rxSymbols] = nrPDSCHDecode(carrier,pdsch,pdschEq,noiseEst);
rxbits = llr{1}<0;

Results

Display the in-phase and quadrature constellation of the received PDSCH symbols and measure their error vector magnitude (EVM). To check for transmission errors, compare the received and transmitted bits.

figure
plot(pdschEq,'o')
xlabel("In-Phase")
ylabel("Quadrature")
title('Equalized PDSCH Constellation')

measureEVM = comm.EVM;
EVM = measureEVM(pdschSymbols,pdschEq);
fprintf("Average PDSCH EVM: %.3g %%.",mean(EVM))

Average PDSCH EVM: 1.3 %.

transmissionError = ~isequal(rxbits,data)

transmissionError = logical
   0

Local Functions

function F  = analogBeamformer(subarraySize,elementSpacing,numRFChains,azimuth,downTilt)
% Define analog beamformers with directions defined by the set of azimuth
% and downtilt input angles.

    % Expand azimuth and elevation as required
    azimuth= azimuth(:);
    downTilt = downTilt(:);
    if isscalar(azimuth)
        azimuth = repmat(azimuth,numel(downTilt),1);
    end    
    if isscalar(downTilt)
        downTilt = repmat(downTilt,numel(azimuth),1);
    end
    
    % Arrange azimuth and downtilt values in the third dimension of the
    % array for convenience
    azimuth = permute(azimuth,[3 2 1]);
    downTilt = permute(downTilt,[3 2 1]);

    % Spherical to Cartesian coordinates factors
    az = azimuth*pi/180;
    el = -downTilt*pi/180;
    r = 1;
    [~,ys,zs] = sph2cart(az,el,r);

    % Vertical steering vector
    Nz = subarraySize(1);
    dz = elementSpacing(1);
    vsv = exp(-1i*2*pi*dz*((0:Nz-1)-(Nz-1)/2).'.*zs);

    % Horizontal steering vector
    Ny = subarraySize(2);
    dy = elementSpacing(2);
    hsv = exp(-1i*2*pi*dy*((0:Ny-1)-(Ny-1)/2).*ys);

    % Subarray steering vector
    F = reshape(pagemtimes(vsv,hsv),[],1);
    
    % Replicate beamformers to match the number of RF chains
    s = min(numel(azimuth),numRFChains);
    F = repmat(F,ceil(numRFChains/s),1);
    F = F(1:numRFChains*Nz*Ny);

end

function RFConnections = RFChainConnections(antArray)
% Define the connections between RF chains and antenna array elements

    % Antenna array dimensions
    M = antArray.Size(1);
    N = antArray.Size(2);
    Mas = antArray.SubarraySize(1);
    Nas = antArray.SubarraySize(2);
    Msa = M/Mas;
    Nsa = N/Nas;
    P = antArray.NumPolarizations;

    if mod(Msa,1) || mod(Nsa,1)
        error("Array and subarray sizes must be such that the array can be partitioned into an integer number of equal-size subarrays.")
    end

    % Determine the antenna element indices connected to each RF chain
    RFConnections = zeros(M,N,P);
    numSubarrays = Msa*Nsa;
    for p = 1:P
        RFConnections(:,:,p) = kron(reshape(1:numSubarrays,Msa,Nsa),ones(Mas,Nas)) + numSubarrays*(p-1);
    end

    % Check that there are no repetitions
    numAntElem = M*N*P;
    for a = 0:numAntElem-1
        rfca = RFConnections(a + (1:numAntElem:end));
        rfca(rfca==0) = NaN;
        if numel(rfca) ~= numel(unique(rfca))
            error("The same RF chain has been connected to the same antenna multiple times.")
        end
    end
end

function outWaveform = analogBeamformWaveform(waveform,F,RFConnections,NumAntElements)
% Apply beamforming weights F to input waveform. The beamforming weights F
% must be sorted following the RF chain connection list RFConnections.

    % Reshape RF connections to 2-D (NumAntElements-by-NRF) if needed. NRF is
    % the number of RF chains connected to an antenna element. For more
    % information, see the RF Connections to Antenna Elements section.
    if ~(height(RFConnections) == NumAntElements) || ~ismatrix(RFConnections)
        RFConnections = reshape(RFConnections,NumAntElements,[]);
    end
    
    % Preallocate output waveform (NTimeSamples-by-NumAntElements)
    NTimeSamples = size(waveform,1);
    outWaveform = zeros(NTimeSamples,NumAntElements);

    % Apply beamforming weights to input waveform
    NRF = size(waveform,2);
    for rf = 1:NRF
        rfAntennaElement = RFConnections == rf;
        antennaElement = any(rfAntennaElement,2);
        outWaveform(:,antennaElement) = outWaveform(:,antennaElement) + waveform(:,rf).*F(rfAntennaElement).' ;
    end

end

function estChannelGrid = combineChannelEstimate(estChannelGrid,F,RFConnections,NumAntElements)
% Combine input perfect channel estimates with input analog beamforming
% weights F. The beamforming weights F must be sorted following the RF
% chain connections RFConnections. The input perfect channel estimate array
% of size [K-by-N-by-Nr-by-Nt] contains estimates for each transmit antenna
% element. The output array of size [K-by-N-by-Nr-by-NRF] contains perfect
% channel estimates for each RF chain.

    % Reshape RF connection list to 2-D (NumAntennaElements-by-N) if needed
    if ~(height(RFConnections) == NumAntElements) || ~ismatrix(RFConnections)
        RFConnections = reshape(RFConnections,NumAntElements,[]);
    end

    NRF = length(unique(RFConnections));

    % Change estimate dimension order for convenience before combining
    H = permute(estChannelGrid,[3 4 1 2]);

    % Combine estimates using beamforming weights
    Hbf = zeros([size(H,1),NRF,size(H,[3,4])]);
    for rf = 1:NRF
        rfAntennaElement = RFConnections == rf;
        antennaElement = any(rfAntennaElement,2);
        Hbf(:,rf,:,:) = pagemtimes(H(:,antennaElement,:,:),F(rfAntennaElement));
    end

    % Change estimate dimension order back
    estChannelGrid = permute(Hbf,[3 4 1 2]);
end

function channel = createChannel(simParameters)
% Create and configure the propagation channel    
    
    if contains(simParameters.DelayProfile,'CDL')

        % Create CDL channel
        channel = nrCDLChannel;

        % Tx antenna array configuration in CDL channel. The size of the
        % antenna array is [M,N,P,Mg,Ng]. M and N are the number of rows
        % and columns in the antenna array, respectively. P is the number
        % of polarizations (1 or 2). Mg and Ng are the number of row and
        % column array panels, respectively.
        txArray = simParameters.TransmitAntennaArray; 

        channel.TransmitAntennaArray.Size = [txArray.Size txArray.NumPolarizations 1 1];
        channel.TransmitAntennaArray.ElementSpacing = [txArray.ElementSpacing 1 1]; % Element spacing in wavelengths
        channel.TransmitAntennaArray.PolarizationAngles = [-45 45];                 % Polarization angles in degrees
        
        % Rx antenna array configuration in CDL channel
        rxArray = simParameters.ReceiveAntennaArray;

        channel.ReceiveAntennaArray.Size = [rxArray.Size rxArray.NumPolarizations 1 1];
        channel.ReceiveAntennaArray.ElementSpacing = [rxArray.ElementSpacing 1 1];  % Element spacing in wavelengths
        channel.ReceiveAntennaArray.PolarizationAngles = [0 90];                    % Polarization angles in degrees

    else
        error('Channel not supported.')
    end

    % Configure other channel parameters: delay profile, delay spread, and
    % maximum Doppler shift
    channel.DelayProfile = simParameters.DelayProfile;
    channel.DelaySpread = simParameters.DelaySpread;
    channel.MaximumDopplerShift = simParameters.MaximumDopplerShift;

    % Get information about the baseband waveform after OFDM modulation step
    waveInfo = nrOFDMInfo(simParameters.Carrier);

    % Update channel sample rate based on carrier information
    channel.SampleRate = waveInfo.SampleRate;

    % Specify the channel response output to obtain the OFDM response of
    % the channel
    channel.ChannelResponseOutput = "ofdm-response";
    
end

function [channel,maxChannelDelay] = setupChannel(simParameters)
% Reset the propagation channel and obtain the maximum channel delay

    % Extract channel
    channel = simParameters.Channel;
    channel.reset();

    % Get the channel information
    chInfo = info(channel);
    maxChannelDelay = chInfo.MaximumChannelDelay;

end

function [N0, noiseEst] = setupReceiver(simParameters,channel)
% Calculate noise standard deviation and noise variance

    % Calculate noise standard deviation. Normalize noise power by the FFT
    % size used in OFDM modulation, as the OFDM modulator applies this
    % normalization to the transmitted waveform.
    SNRdB = simParameters.SNRIn;
    SNR = 10^(SNRdB/10);
    waveInfo = nrOFDMInfo(simParameters.Carrier);
    N0 = 1/sqrt(double(waveInfo.Nfft)*SNR);
    
    % Also normalize by the number of receive antennas if the channel
    % applies this normalization to the output
    if channel.NormalizeChannelOutputs
        chInfo = info(channel);
        N0 = N0/sqrt(chInfo.NumReceiveAntennas);
    end

    noiseEst = N0^2*double(waveInfo.Nfft);
end

function numElemenets = numAntennaElements(antArray)
% Calculate number of antenna elements in an antenna array

    numElemenets = antArray.NumPolarizations*prod(antArray.Size);
    
end

function validateNumLayers(simParameters)
% Validate the number of layers, relative to the antenna geometry

    numlayers = simParameters.PDSCH.NumLayers;
    numRFChains = prod(simParameters.TransmitAntennaArray.NumRFChains);
    nrxants = simParameters.NRxAnts;
    antennaDescription = sprintf('min(numRFChains,NRxAnts) = min(%d,%d) = %d',numRFChains,nrxants,min(numRFChains,nrxants));
    if numlayers > min(numRFChains,nrxants)
        error('The number of layers (%d) must satisfy NumLayers <= %s', ...
            numlayers,antennaDescription);
    end

    % Display a warning if the maximum possible rank of the channel equals
    % the number of layers
    if (numlayers > 2) && (numlayers == min(numRFChains,nrxants))
        warning(['The maximum possible rank of the channel, given by %s, is equal to NumLayers (%d).' ...
            ' This can result in a decoding failure in certain channel conditions.' ...
            ' Try decreasing the number of layers or increasing the channel rank' ...
            ' (use more transmit or receive antennas).'],antennaDescription,numlayers); %#ok<SPWRN>
    end

end

function plotRFConnectionsToSubarrays(M,antArray)
% Plot RF connections to subarrays

    clims = [min(M(:)) max(M(:))];
    if clims(1) == clims(2)
        clims(2) = clims(1)+1;
    end

    for i = 1:size(M,3)
        figure
        imagesc(M(:,:,i),clims); 

        % Plot lines separating elements and subarrays
        W = size(M,2);
        H = size(M,1);
        Ws = antArray.SubarraySize(2);
        Hs = antArray.SubarraySize(1);
        plotLines(W,H,Ws,Hs,3,'w')
        plotLines(W,H,1,1,0.5,'w')
        
        % Quantize colormap
        cm = colormap;
        colormap(cm(1:floor(size(cm,1)/(clims(2)-clims(1))-1):end,:)); % Adjust colormap to discrete values
        cb = colorbar;
        ylabel(cb,'RF Chain');

        % Add title and labels
        title("RF Connections to Antenna Elements (Pol = " + i + ")");
        xlabel('Antenna Element (H)')
        ylabel('Antenna Element (V)')
        set(gca,XTick = 1:W)
        set(gca,YTick = 1:H)
        hold on;
    end
    
end

function plotLines(W,H,Ws,Hs,lw,col)

    for row = 1:Hs:H+1
        line(0.5+[0,W],-0.5+[row,row],'Color',col,'LineWidth',lw);
    end
    for column = 1:Ws:W+1
        line([column,column]-0.5,0.5+[0,H],'Color',col,'LineWidth',lw);
    end

end

function plotBeams(subarrayDims,elementSpacing,F,rfIndex)

    if nargin == 3
        rfIndex = 1;
    end
    
    % Define 3-D plotting space in Cartesian coordinates
    ph = linspace(0,2*pi,100);
    th = linspace(-pi/2,pi/2,120);

    [TH,PH] = ndgrid(th,ph);
    [x,y,z] = sph2cart(PH,TH,1);
    r = [x(:) y(:) z(:)];
    
    % Calculate positions of each antenna element in a subarray
    nh = subarrayDims(2);
    dh = elementSpacing(2);
    hvec = dh*((0:nh-1)-(nh-1)/2);

    nv = subarrayDims(1);
    dv = elementSpacing(1);
    vvec = dv*((0:nv-1)-(nv-1)/2);

    [Dx,Dy] = meshgrid(hvec,vvec);
    d = [zeros(numel(Dx),1) Dx(:) Dy(:)];
    
    % Plot pattern using the input beam weights
    figure
    for rf = 1:length(rfIndex)
        thisF = F((rfIndex(rf)-1)*nh*nv + (1:nh*nv));
        A = reshape(sum(thisF.*exp(1i*2*pi*d*r'),1), size(x));
        [Ax,Ay,Az] = sph2cart(PH,TH,abs(A));
        nexttile;
        surf(Ax,Ay,Az,abs(A),'EdgeAlpha',0.1); axis equal
        title("Subarray Pattern " + rf);
        xlabel('x');
        ylabel('y');
    end
    colorbar

end