5G 波形を使用したセンシングと通信の統合
この例では、TR 38.901 ISAC チャネル モデルを使用して、5G New Radio (NR) 波形のセンシングと通信の統合 (ISAC) 機能について説明します。物理ダウンリンク共有チャネル (PDSCH) を通じた 5G NR リンクのデータ フレーム伝送をシミュレーションします。また、受信フレームから取得されたチャネル行列推定値を処理し、チャネル内に存在する移動中の無人航空機 (UAV) のセンシング ターゲット (ST) のレーダー測定を抽出する方法を示します。
はじめに
ISAC の中核となる概念は、同一のハードウェアと周波数スペクトルを利用して、通信機能とセンシング機能の両方を実行することです。近年の研究では、新たなデュアル機能システムの共同設計から、既存の無線ネットワークでのセンシング機能の実現に至るまで、さまざまな ISAC アプローチが検討されています。後者のアプローチは、既存の装置やネットワーク インフラストラクチャを活用するため、特に有効です。Beyond 5G (B5G) や 6G などの次世代無線ネットワークにセンシングを組み込むことで、電磁スペクトルの混雑を緩和できるだけでなく、各種の新たな用途を実現できます。これには、高精度な位置特定やトラッキング、さらには人の行動認識が含まれます。
この例では、3GPP NR 規格で定義された 5G NR 波形を使用して、移動中の UAV センシング ターゲットのセンシングを効果的に実現する方法を示します。受信した PDSCH フレームから取得されるチャネル行列推定値には、送信波形が受信機に到達する過程で生じる時間遅延とドップラー シフトの情報が元々含まれています。この情報を標準的なレーダー データ キューブ手法と同様の方法で処理することで、これらの時間遅延とドップラー シフトを対応する UAV センシング ターゲットの位置と速度に変換できます。これによりレーダー検出が形成されます。
この例では、単一のセンシング送信機 (STX) と単一のセンシング受信機 (SRX) を使った ISAC-UAV センシング シナリオを定義することから始めます。TRP-UE バイスタティック センシング モードを適用します。このモードで、STX は基地局を表し、SRX はユーザー端末をモデル化します。この例では、次に、送信する 5G NR 波形を構成し、復調基準信号 (DM-RS) のパラメーターの選択が目的とするセンシング パフォーマンスと直接関係していることを示します。続いて、この例では、TR 38.901 ISAC チャネルを通じた PDSCH フレームの伝搬をシミュレーションします。このチャネルは、センシング ターゲットと背景環境の両方の影響を受けるマルチパス成分をモデル化します。各受信フレームについて、チャネル行列を推定し、その後、チャネル内に存在する移動中の UAV ST を検出するために処理します。
この例では、センシングの観点から、モデル化された 5G リンクが、ターゲットのバイスタティック範囲とバイスタティック到来角を測定できるバイスタティック レーダーであることを示します。最後に、ターゲットのトラックを形成し、直交座標におけるターゲットの位置と速度を推定するために、これらの測定値がトラッキング アルゴリズムに渡されます。
ISAC シナリオ
周波数 30 GHz で動作する 5G NR リンクについて考えます。
rng('default'); % Set the random number generator for reproducibility carrierFrequency = 30e9; % Carrier frequency (Hz) wavelength = freq2wavelen(carrierFrequency);
STX (gNodeB) が xy 平面の原点にあり、高さが 10 メートルであるものとします。
txPosition = [0; 0; 10]; % Base station locationSTX には、半波長間隔で配置された 8 つの等方性アンテナ素子から成る等間隔直線アレイ (ULA) があるものとします。送信機電力を 46 dBm に設定します。
numTxAntennas = 8; % Number of antenna elements at the Tx array element = phased.IsotropicAntennaElement('BackBaffled', true); txArray = phased.ULA(numTxAntennas, wavelength/2, 'Element', element); txPower = 46; % dBm
送信アレイの法線を x 軸方向に向けます。txOrientation ベクトルは、TR 38.901 の Section 7.1.3 で規定されているように、ベアリング、ダウンチルト、およびスラントの回転角度を度単位でそれぞれ指定します。
txOrientation = [0; 0; 0]; % Tx array orientation [bearing; downtilt; slant] (degrees)SRX (ユーザー端末) が送信機からある程度離れた位置にあるものとします。
rxPosition = [50; 60; 2]; % User equipment locationSRX には、半波長間隔で配置された 8 つの等方性アンテナ素子から成る等間隔直線アレイ (ULA) があるものとします。受信機のノイズ指数を 5 dB に設定します。
numRxAntennas = 8; % Number of antenna elements at the Rx array rxArray = phased.ULA(numRxAntennas, wavelength/2, 'Element', element); rxNoiseFigure = 5; % dB
受信アレイの法線を負の y 軸方向に向けます。
rxOrientation = [270; 0; 0]; % Rx array orientation [bearing; downtilt; slant] (degrees)helperGetTargetTrajectories 補助関数を使用し、移動中の 2 つの UAV センシング ターゲットに対して、直線の等速軌道を 2 つ生成します。これらは、センシング機能の対象ターゲットです。
targets.Trajectories = helperGetTargetTrajectories();
ans=2×5 table
Trajectory Start Position [x, y, z] End Position [x, y, z] Length (m) Speed (m/s)
__________ ________________________ ______________________ __________ ___________
"1" {[25 -5 10]} {[40 30 10]} 38.079 9.5197
"2" {[70 30 10]} {[80 10 10]} 22.361 5.5902
numTargets = numel(targets.Trajectories);
h38901ISACChannel System object (TM) オブジェクトを使用して、STX、ST、および SRX の間のチャネルをモデル化します。
channel = h38901ISACChannel; channel.SensingScenario = "ISAC-UAV"; channel.SensingMode = "TRP-UE Bistatic"; channel.CommunicationScenario = "UMi"; % Sensing Transmitter channel.STX.Position = txPosition; channel.STX.TransmitAntennaArray = txArray; channel.STX.TransmitArrayOrientation = txOrientation; % Sensing Receiver channel.SRX.Position = rxPosition; channel.SRX.ReceiveAntennaArray = rxArray; channel.SRX.ReceiveArrayOrientation = rxOrientation; % Sensing Targets channel.STs = repmat(channel.STs,1,numTargets); channel.STs(1) = applyInitialPose(channel.STs(1),targets.Trajectories{1}); channel.STs(2) = applyInitialPose(channel.STs(2),targets.Trajectories{2}); % Configure channel properties channel.CenterFrequency = carrierFrequency; channel.LOSProbability = 1; channel.TargetPowerOffset = 15; % (dB) for power normalization channel.Seed = 3; channel.ScenarioExtents = helperGetScenarioExtents(channel,targets.Trajectories);
helperVisualizeScenario 補助関数を使用して、シナリオを可視化します。
helperVisualizeScenario(channel,targets.Trajectories);
title('ISAC Scenario');
センシングの観点では、このシナリオにおいて送信機と受信機がバイスタティック ペアを形成することに注意してください。
5G 波形の構成
5G NR のシステム パラメーターを設定します。チャネル帯域幅を 50 MHz とします。
channelBandwidth = 50; % Transmission bandwidth and guardbands (MHz)帯域幅占有率を設定します。これは、伝送帯域幅とチャネル帯域幅の比率です。帯域幅占有率を低くすると、保護帯域のサイズが大きくなります。これにより、スペクトル効率が低下しますが、スペクトル放射は低減されます。
bandwidthOccupancy = 0.9; % Ratio of transmission bandwidth to channel bandwidthサブキャリア間隔を 60 kHz とします。
subcarrierSpacing = 60; % Subcarrier spacing (kHz)helperGet5GWaveformConfiguration を使用して 5G 波形の構成を作成します。
waveformConfig = helperGet5GWaveformConfiguration(channelBandwidth, bandwidthOccupancy, subcarrierSpacing); waveformInfo = nrOFDMInfo(waveformConfig.Carrier)
waveformInfo = struct with fields:
Nfft: 1024
SampleRate: 61440000
CyclicPrefixLengths: [104 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 104 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72 72]
SymbolLengths: [1128 1096 1096 1096 1096 1096 1096 1096 1096 1096 1096 1096 1096 1096 1096 1096 1096 1096 1096 1096 1096 1096 1096 1096 1096 1096 1096 1096 1128 1096 1096 1096 1096 1096 1096 1096 1096 1096 1096 1096 1096 1096 … ] (1×56 double)
Windowing: 36
SymbolPhases: [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
SymbolsPerSlot: 14
SlotsPerSubframe: 4
SlotsPerFrame: 40
transmissionBandwidth = subcarrierSpacing*1e3*waveformConfig.Carrier.NSizeGrid*12; % Bandwidth of resource blocks (Hz) fprintf("Transmission bandwidth: %.2f MHz\n",transmissionBandwidth*1e-6);
Transmission bandwidth: 44.64 MHz
チャネルと受信機のサンプル レートを波形のサンプル レートと同じに設定します。
channel.SampleRate = waveformInfo.SampleRate; receiver.SampleRate = waveformInfo.SampleRate;
PDSCH フレームには、受信機がチャネル特性を推定するのに役立つ DM-RS 信号が含まれています。通信機能は、この情報を使用して送信データを正しく復調します。センシング機能は、同じ情報を使用してチャネル内に存在するセンシング ターゲットの位置と速度を推定できます。DM-RS パイロット信号の構成は、センシング パフォーマンスに直接影響します。DM-RS のさまざまな構成の詳細については、NR PDSCH のリソース割り当てと DM-RS および PT-RS 基準信号 (5G Toolbox)の例を参照してください。
DM-RS のパラメーターを設定します。
waveformConfig.PDSCH.DMRS.DMRSTypeAPosition = 2; % Mapping type A only. First DM-RS symbol position is 2 waveformConfig.PDSCH.DMRS.DMRSLength = 1; % Single-symbol DM-RS waveformConfig.PDSCH.DMRS.DMRSAdditionalPosition = 2; % Additional DM-RS symbol positions (max range 0...3) waveformConfig.PDSCH.DMRS.NumCDMGroupsWithoutData = 1; % Number of CDM groups without data waveformConfig.PDSCH.DMRS.DMRSConfigurationType = 1; % DM-RS configuration type (1,2)
構成済みリソース グリッドのうちの 1 つのリソース ブロックを可視化します。
resourceGrid = nrResourceGrid(waveformConfig.Carrier, waveformConfig.PDSCH.NumLayers);
ind = nrPDSCHDMRSIndices(waveformConfig.Carrier, waveformConfig.PDSCH);
resourceGrid(ind) = nrPDSCHDMRS(waveformConfig.Carrier, waveformConfig.PDSCH);
helperVisualizeResourceGrid(abs(resourceGrid(1:12,1:14,1)));
title({'PDSCH DM-RS Resource Elements in the Carrier Resource Grid', '(single resource block)'});
時間領域におけるパイロット信号の密度は、エイリアシングを生じさせることなく、明瞭に測定できる最大ドップラー シフトを決定します。構成済みの PDSCH フレームでは、時間領域におけるパイロット信号間の最大間隔が = 5 シンボルとなります。対応する最大明瞭ドップラー シフトは、次の式を使用して計算できます。
ここで、 はシンボル区間です。
% Compute maximum pilot signal separation in time [~,pdschInfo] = nrPDSCHIndices(waveformConfig.Carrier,waveformConfig.PDSCH); Mt = max(diff(pdschInfo.DMRSSymbolSet)); % Duration of an OFDM symbol in seconds ofdmSymbolDuration = max(waveformInfo.SymbolLengths)/waveformInfo.SampleRate; maxDoppler = 1/(2*Mt*ofdmSymbolDuration); fprintf('Maximum unambiguous Doppler shift: %.2f kHz\n', maxDoppler*1e-3);
Maximum unambiguous Doppler shift: 5.45 kHz
対応する最大速度を計算し、検討対象の ISAC シナリオにおけるターゲットの最大対地速度と比較します。
maxVelocity = dop2speed(maxDoppler, wavelength)/2;
fprintf('Maximum unambiguous velocity: %.2f m/s\n', maxVelocity);Maximum unambiguous velocity: 27.22 m/s
このシナリオにおけるターゲットの速度は、最大明瞭速度の範囲内にあります。
同様に、周波数領域におけるパイロット信号の密度は、エイリアシングを生じさせることなく測定できる最大時間遅延を決定します。周波数領域では、パイロット信号が = 2 サブキャリアごとに繰り返されます。対応する最大明瞭時間遅延は、次の式を使用して計算できます。
ここで、 はサブキャリア間隔です。
% Compute pilot signal separation in frequency Mf = min(diff(waveformConfig.PDSCH.DMRS.DMRSSubcarrierLocations)); maxDelay = 1/(2*Mf*subcarrierSpacing*1e3); fprintf('Maximum time delay: %.2fe-6 s\n', maxDelay*1e6);
Maximum time delay: 4.17e-6 s
対応する最大レンジを計算します。
% Baseline distance between Tx and Rx baseline = vecnorm(txPosition - rxPosition); % Propagation speed propagationSpeed = physconst("LightSpeed"); % Bistatic maximum range maxRange = maxDelay*propagationSpeed-baseline; fprintf('Maximum unambiguous bistatic range: %.2f m\n', maxRange);
Maximum unambiguous bistatic range: 1170.62 m
このシナリオのすべてのターゲットは、最大明瞭レンジの十分な範囲内にあります。
最後に、伝送帯域幅はレンジ分解能を決定します。
% Range resolution rangeResolution = propagationSpeed/transmissionBandwidth; fprintf('Range resolution: %.2f m\n', rangeResolution);
Range resolution: 6.72 m
システムのレンジ分解能未満のレンジでしか離れておらず、同じ到来角を共有するターゲットは、区別できません。検討対象のシナリオでは、対象となる 2 つのターゲットがレンジと角度の両方で十分に離れています。
チャネル シミュレーションとセンシング データ処理
この例では、10 個の PDSCH フレームの送信をシミュレーションします。各フレームは 0.4 秒間隔で送信され、合計 3.6 秒の時間範囲をカバーします。この手法は、ターゲットの位置の変化を効果的に示すために必要です。この時間間隔内のすべてのフレームをシミュレーションすると、計算量が非常に大きくなります。ただし、シミュレーションされたフレームの間も、5G NR リンクが通常どおり動作し続けるものと仮定します。
% Number of PDSCH frames to simulate numSensingFrames = 10; % Time step between simulated PDSCH frames (s) dt = 0.4; % Simulation times t = (0:numSensingFrames-1)*dt;
チャネル行列推定値は、各 PDSCH フレーム伝送から取得します。チャネル行列の次元は です。ここで、 はサブキャリアの数、 は 1 フレーム内の OFDM シンボルの総数、 は受信アンテナの数、 は基準信号ポートの数です。この例では、基地局はビームフォーミングを実行しません。そのため、 は送信アンテナの数と等しくなります。チャネル行列推定値はレーダー データ キューブと見なすことができます。各次元は次を表します。
サブキャリア次元はレンジ (ファストタイム)
シンボル次元はスロータイム
受信アンテナ次元は到来角 (AoA)
基準信号ポート次元は発射角 (AoD)
チャネル行列推定値の各次元は個別に処理できます。その後、結果はスロータイム次元と AoD 次元にわたって非コヒーレントに積分され、レンジ-AoA マップが生成されます。これはターゲットの位置を推定するために使用されます。
この処理を実行するには、AoA のグリッドを定義し、対応するステアリング ベクトルを計算します。
numAoA = ceil(2*pi*numRxAntennas); % AoA grid size aoaGrid = linspace(-90, 90, numAoA); % AoA grid % Receive array steering vectors rxSteeringVector = phased.SteeringVector(SensorArray=channel.SRX.ReceiveAntennaArray); rxsv = rxSteeringVector(carrierFrequency, aoaGrid);
AoD のグリッドを定義し、対応するステアリング ベクトルを計算します。
numAoD = ceil(2*pi*numTxAntennas); % AoD grid size aodGrid = linspace(-90, 90, numAoD); % AoD grid % Transmit array steering vectors txSteeringVector = phased.SteeringVector(SensorArray=channel.STX.TransmitAntennaArray); txsv = txSteeringVector(carrierFrequency, aodGrid);
バイスタティック範囲グリッドと対応するファストタイム ステアリング ベクトルを定義します。
% Limit the maximum bistatic range of interest to 100 meters rangeGrid = -20:rangeResolution/4:100; % Bistatic range grid numRange = numel(rangeGrid); % Bistatic range grid size % Fast-time steering vectors rangeFFTBins = (rangeGrid+baseline)/rangeResolution; numSubcarriers = waveformConfig.Carrier.NSizeGrid*12; k = (0:numSubcarriers-1).'; rngsv = (exp(2*pi*1i*k*rangeFFTBins/numSubcarriers)/numSubcarriers);
バイスタティック範囲-AoA マップを処理する 2 次元の定誤警報率 (CFAR) 検出器を設定します。レンジ分解能および角度分解能に基づいて、CFAR ガード バンド サイズを決定します。
% Number of guard cells in range numGuardCellsRange = ceil(rangeResolution/(rangeGrid(2)-rangeGrid(1))) + 1; % Resolution in AoA domain aoaResolution = ap2beamwidth((numRxAntennas-1)*channel.SRX.ReceiveAntennaArray.ElementSpacing, wavelength); % Number of guard cells in AoA numGuardCellsAoA = ceil(aoaResolution/(aodGrid(2)-aodGrid(1))) + 1; % Create a 2-D CFAR detector cfar2D = phased.CFARDetector2D('GuardBandSize',[numGuardCellsRange numGuardCellsAoA], 'TrainingBandSize', [4 6],... 'ProbabilityFalseAlarm', 1e-3, 'OutputFormat', 'CUT result', 'ThresholdOutputPort', false); % Compute indices of the cells under test offsetIdxs = cfar2D.TrainingBandSize + cfar2D.GuardBandSize; rangeIdxs = offsetIdxs(1)+1:numel(rangeGrid)-offsetIdxs(1); angleIdxs = offsetIdxs(2)+1:numel(aoaGrid)-offsetIdxs(2); [sumRangeIdxs, aoaIdxs] = meshgrid(rangeIdxs, angleIdxs); cutidx = [sumRangeIdxs(:).'; aoaIdxs(:).'];
CFAR アルゴリズムは、真のターゲット位置の周囲で複数の検出を生成する場合があります。通常、これらの検出はまとめてクラスタリングされ、単一の検出が得られます。DBSCAN クラスタラー オブジェクトを作成します。
clusterer = clusterDBSCAN();
helperConfigureTracker 補助関数を使用してターゲット トラッキング オブジェクトを作成します。返されるトラッカーは、バイスタティック範囲とバイスタティック到来角の測定値を使用してトラッキングを実行するように構成されています。トラッカーの状態ベクトルは直交座標です。したがって、トラッカーはターゲットの位置推定を向上させるだけでなく、ターゲットの速度も推定します。
% Define transmitter and receiver orientation axes txOrientationAxes = rotationMatrix(txOrientation); rxOrientationAxes = rotationMatrix(rxOrientation); % Configure tracker tracker = helperConfigureTracker(txPosition, rxPosition, txOrientationAxes, rxOrientationAxes, rangeResolution, aoaResolution); % Preallocate space to store estimated target state vectors [x; vx; y; vy] % Preallocate space for 10 tracks targetStateEstimates = zeros(4, numSensingFrames, 10); trackIDs = [];
ISAC シナリオ、真のターゲット位置、ターゲット検出、およびトラックをまとめて表示する可視化を生成するには、generateVisualizations を true に設定します。
% Produce visualizations generateVisualizations =false; if generateVisualizations scenarioCFARResultsFigure = figure; tl = tiledlayout(scenarioCFARResultsFigure, 1, 2, 'Padding', 'tight', 'TileSpacing', 'tight');
helperTheaterPlotter を使用して、シナリオ、ターゲット軌跡、検出、およびトラックを可視化します。
scenarioPlotter = helperTheaterPlotter(tl);
scenarioPlotter.plotTxAndRx(txPosition, rxPosition);
scenarioPlotter.plotTrajectories(targets.Trajectories, t);
title(scenarioPlotter, 'ISAC Scenario');helperCFARResultsVisualizer を使用して、CFAR 結果を可視化します。
cfarResultsVisualizer = helperCFARResultsVisualizer(tl);
xlabel(cfarResultsVisualizer, 'AoA (degrees)');
ylabel(cfarResultsVisualizer, 'Bistatic Range (m)');
title(cfarResultsVisualizer, 'CFAR Results');
endシミュレーションを実行し、PDSCH フレームを 1 つずつ送信して処理します。
% Number of symbols per frame numSymbolsPerFrame = waveformInfo.SlotsPerFrame*waveformInfo.SymbolsPerSlot; % Transmit and process one PDSCH frame at a time for i = 1:numSensingFrames fprintf('\nSimulating frame at time %g s\n', t(i)); % helperSimulateLinkSingleFrame returns a 4-D matrix % numSubcarriers x numSymbolsPerFrame x numRxAntennas x numTxAntennas sensingCSI = helperSimulateLinkSingleFrame(channel, waveformConfig, txPower, rxNoiseFigure, targets, t(i)); % Process data cube by putting a zero in the first Doppler bin to filter out background noise X = fft(sensingCSI, [], 2); X(:, 1, :, :) = 0; x = ifft(X, [], 2); % Process data cube along the AoD dimension x = permute(x, [4 1 3 2]); x = txsv'*reshape(x, numTxAntennas, []); % Process data cube along the range dimension x = permute(reshape(x, numAoA, numSubcarriers, numRxAntennas, numSymbolsPerFrame), [2 3 1 4]); x = rngsv.'*reshape(x, numSubcarriers, []); % Process data cube along the AoA dimension x = permute(reshape(x, numRange, numRxAntennas, numAoD, numSymbolsPerFrame), [2 1 3 4]); x = rxsv'*reshape(x, numRxAntennas, []); % Average over AoDs and symbols x = reshape(x, numAoA, numRange, []); sumRangeAoAMap = sqrt(sum(abs(x).^2, 3)).'; % Perform 2-D CFAR detection over range-AoA map cutResult = cfar2D(sumRangeAoAMap, cutidx); detectionIdxs = cutidx(:, cutResult); detectionValues = [rangeGrid(detectionIdxs(1,:)); aoaGrid(detectionIdxs(2,:))]; % Cluster CFAR detections [~, clusterids] = clusterer(detectionValues.'); uniqClusterIds = unique(clusterids); % Loop over detection clusters, and for each cluster find a % center m = numel(uniqClusterIds); clusteredDetections = zeros(2, m); for j = 1:m idxs = clusterids == uniqClusterIds(j); d = detectionValues(:, idxs); % Represent a cluster of detections with its center clusteredDetections(:, j) = mean(d, 2); end % Format detections so that they can be processed by the tracker detections = helperFormatDetectionsForTracker(clusteredDetections, t(i), rangeResolution, aoaResolution); % Pass detections to the tracker tracks = tracker(detections); % Store target state estimates for it = 1:numel(tracks) id = tracks(it).TrackID; trackIDs = union(trackIDs, id); targetStateEstimates(:, i, id) = tracks(it).State; end % Visualization if generateVisualizations % Plot true target positions targetPositions = zeros(numTargets, 3); for it = 1:numTargets targetPositions(it, :) = lookupPose(targets.Trajectories{it}, t(i)); end scenarioPlotter.plotTargetPositions(targetPositions); % Plot target positions estimated from the detections measuredPositions = helperGetCartesianMeasurement(clusteredDetections, txPosition, rxPosition, rxOrientationAxes); scenarioPlotter.plotDetections(measuredPositions); scenarioPlotter.plotTracks(tracks); % Visualize CFAR results [trueRrx, trueAoA] = rangeangle(targetPositions.', rxPosition, rxOrientationAxes); [trueRtx, trueAoD] = rangeangle(targetPositions.', txPosition, txOrientationAxes); trueBistaticRange = trueRrx + trueRtx - baseline; cfarResultsVisualizer.plotCFARImage(aoaGrid, rangeGrid, sumRangeAoAMap); cfarResultsVisualizer.plotTruth([trueBistaticRange; trueAoA(1, :)]); cfarResultsVisualizer.plotClusteredDetections(clusteredDetections); sgtitle(scenarioCFARResultsFigure, sprintf('Frame %d, Simulation Time %.1f s', i, t(i)), 'FontSize' ,10); drawnow; end end
Simulating frame at time 0 s
Throughput(Mbps) for 1 frame(s) = 1.7928 Throughput(%) for 1 frame(s) = 100.0000
Simulating frame at time 0.4 s
Throughput(Mbps) for 1 frame(s) = 1.7928 Throughput(%) for 1 frame(s) = 100.0000
Simulating frame at time 0.8 s
Throughput(Mbps) for 1 frame(s) = 1.7928 Throughput(%) for 1 frame(s) = 100.0000
Simulating frame at time 1.2 s
Throughput(Mbps) for 1 frame(s) = 1.7928 Throughput(%) for 1 frame(s) = 100.0000
Simulating frame at time 1.6 s
Throughput(Mbps) for 1 frame(s) = 1.7928 Throughput(%) for 1 frame(s) = 100.0000
Simulating frame at time 2 s
Throughput(Mbps) for 1 frame(s) = 1.7928 Throughput(%) for 1 frame(s) = 100.0000
Simulating frame at time 2.4 s
Throughput(Mbps) for 1 frame(s) = 1.7928 Throughput(%) for 1 frame(s) = 100.0000
Simulating frame at time 2.8 s
Throughput(Mbps) for 1 frame(s) = 1.7928 Throughput(%) for 1 frame(s) = 100.0000
Simulating frame at time 3.2 s
Throughput(Mbps) for 1 frame(s) = 1.7928 Throughput(%) for 1 frame(s) = 100.0000
Simulating frame at time 3.6 s
Throughput(Mbps) for 1 frame(s) = 1.7928 Throughput(%) for 1 frame(s) = 100.0000
シミュレーション結果は、チャネル行列推定値を処理することで、バイスタティック範囲-AoA 空間でレーダー検出を正常に形成できることを示しています。これらの検出は、ターゲットの真の位置に近接しています。その後、トラッカーはこれらの検出を利用して、ターゲットの状態を推定します。これは、直交座標における位置と速度で構成されます。
helperPlotSpeedEstimationResults 補助関数を使用して、トラッカー出力に基づくターゲットの速度と進行方向の推定値を真の値と比較します。
helperPlotSpeedEstimationResults(t, targets.Trajectories, targetStateEstimates, trackIDs)
トラックを開始し、速度の計算を始めるには、トラッカーに少なくとも 2 つの初期検出が必要であることに注意してください。そのため、速度と進行方向の推定値が真の値に収束するまでには、複数回のトラック更新が必要になる場合があります。
まとめ
この例では、PDSCH フレームの送信中に取得されるチャネル行列推定値からバイスタティック レーダー情報を抽出するプロセスを示すことで、5G NR 波形のセンシング機能について説明します。選択した DM-RS パラメーターがセンシング パフォーマンスに与える影響を考察します。さらに、この例では、複数の PDSCH フレームを送信する ISAC シナリオをシミュレーションします。各受信フレームについて、チャネル行列はレーダー データ キューブと同様の方法で処理されます。例では、バイスタティック範囲とバイスタティック AoA の推定のための信号処理チェーンの設定方法を示します。この処理チェーンには、レンジ応答と角度応答の計算に加え、それに続く CFAR 検出と DBSCAN クラスタリングが含まれており、これらにより検出が生成されます。最後に、バイスタティック範囲-AoA 空間における検出がトラッキング アルゴリズムに入力され、ターゲットの位置と速度の推定値が直交座標で計算されます。
サポート関数
function tracker = helperConfigureTracker(txPosition, rxPosition, txOrientationAxes, rxOrientationAxes, rangeResolution, aoaResolution) % Sensor specification sensorSpec = trackerSensorSpec("aerospace", "radar", "bistatic"); sensorSpec.MeasurementMode = "range-angle"; sensorSpec.IsReceiverStationary = true; sensorSpec.IsEmitterStationary = true; sensorSpec.HasElevation = false; sensorSpec.HasRangeRate = false; sensorSpec.MaxNumLooksPerUpdate = 1; sensorSpec.MaxNumMeasurementsPerUpdate = 10; sensorSpec.EmitterPlatformPosition = txPosition; sensorSpec.EmitterPlatformOrientation = txOrientationAxes.'; sensorSpec.ReceiverPlatformPosition = rxPosition; sensorSpec.ReceiverPlatformOrientation = rxOrientationAxes.'; sensorSpec.RangeResolution = rangeResolution; sensorSpec.AzimuthResolution = aoaResolution; sensorSpec.ReceiverFieldOfView = [360 180]; sensorSpec.EmitterFieldOfView = [360 180]; sensorSpec.ReceiverRangeLimits = [0 500]; sensorSpec.EmitterRangeLimits = [0 500]; sensorSpec.DetectionProbability = 0.9; % Target specification targetSpec = trackerTargetSpec('custom'); targetSpec.StateTransitionModel = targetStateTransitionModel('constant-velocity'); targetSpec.StateTransitionModel.NumMotionDimensions = 2; targetSpec.StateTransitionModel.VelocityVariance = 15^2/3*eye(2); targetSpec.StateTransitionModel.AccelerationVariance = 0.01^2/3*eye(2); tracker = multiSensorTargetTracker(targetSpec,sensorSpec,"jipda"); tracker.ConfirmationExistenceProbability = 0.98; tracker.MaxMahalanobisDistance = 10; end function detections = helperFormatDetectionsForTracker(clusteredDetections, currentTime, rangeResoultion, aoaResolution) numDetections = size(clusteredDetections, 2); detections.ReceiverLookTime = currentTime; detections.ReceiverLookAzimuth = 0; detections.ReceiverLookElevation = 0; detections.EmitterLookAzimuth = 0; detections.EmitterLookElevation = 0; detections.DetectionTime = ones(1, numDetections) * currentTime; detections.Azimuth = clusteredDetections(2, :); detections.Range = clusteredDetections(1, :); % The measurement accuracy is proportional to the sensor resolution and % depends on the SNR. detections.AzimuthAccuracy = ones(1, numDetections) * aoaResolution/12; detections.RangeAccuracy = ones(1, numDetections) * rangeResoultion/8; end function meascart = helperGetCartesianMeasurement(dets, txPosition, rxPosition, rxOrientationAxes) if ~isempty(dets) n = size(dets, 2); meassph = zeros(3, n); meassph(1, :) = dets(2, :); meassph(3, :) = dets(1, :); meassph = local2globalcoord(meassph, 'ss', [0; 0; 0], rxOrientationAxes); meascart = bistaticposest(meassph(3, :), meassph(1:2, :), eps*ones(1, n), repmat([eps; eps], 1, n),... txPosition, rxPosition, 'RangeMeasurement', 'BistaticRange'); else meascart = []; end end function trajectories = helperGetTargetTrajectories() % Target 1 start and end points (x,y,z) waypoints = [25 -5 10; 40 30 10]; timeOfArrival = [0 4]; traj1 = waypointTrajectory(waypoints,timeOfArrival); % Target 2 start and end points (x,y,z) waypoints = [70 30 10; 80 10 10]; timeOfArrival = [0 4]; traj2 = waypointTrajectory(waypoints,timeOfArrival); trajectories = {traj1, traj2}; Trajectory = strings(numel(trajectories), 1); Speed = zeros(numel(trajectories), 1); Length = zeros(numel(trajectories), 1); Start = cell(numel(trajectories), 1); Stop = cell(numel(trajectories), 1); for i = 1:numel(trajectories) Trajectory(i) = i; Start{i} = trajectories{i}.Waypoints(1,:); Stop{i} = trajectories{i}.Waypoints(end,:); Length(i) = vecnorm(trajectories{i}.Waypoints(1, :) - trajectories{i}.Waypoints(end, :)); Speed(i) = max(trajectories{i}.GroundSpeed); end table(Trajectory, Start, Stop, Length, Speed, 'VariableNames',... {'Trajectory', 'Start Position [x, y, z]', 'End Position [x, y, z]', 'Length (m)', 'Speed (m/s)'}) end function extents = helperGetScenarioExtents(channel,trajectories) % Derive channel's scenario extents from target trajectories and STX/SRX positions wp = cell(1,numel(trajectories)); for i = 1:numel(trajectories) wp{i} = trajectories{i}.Waypoints; end wp = cat(3,wp{:}); txPos = channel.STX.Position; rxPos = channel.SRX.Position; X = [squeeze(wp(1,1,:));squeeze(wp(2,1,:));txPos(1);rxPos(1)]; Y = [squeeze(wp(1,2,:));squeeze(wp(2,2,:));txPos(2);rxPos(2)]; extentMargin = 1; % m minX = min(X) - extentMargin; minY = min(Y) - extentMargin; width = max(X) - minX + extentMargin; height = max(Y) - minY + extentMargin; extents = [minX minY width height]; end function helperPlotSpeedEstimationResults(t, trajectories, targetStateEstimates, trackIDs) fig = figure; tl = tiledlayout(fig, 1, 2); speedAxes = nexttile(tl); hold(speedAxes, 'on'); xlabel(speedAxes, 'Time (s)'); ylabel(speedAxes, 'Speed (m/s)'); grid(speedAxes, 'on'); title(speedAxes, 'Speed'); xlim(speedAxes, [t(1) t(end)]); lgd = legend(speedAxes, 'Location', 'southeast', 'Orientation', 'horizontal', 'NumColumns', 2); lgd.Layout.Tile = 'south'; headingAxes = nexttile(tl); hold(headingAxes, 'on'); xlabel(headingAxes, 'Time (s)'); ylabel(headingAxes, 'Heading (deg)'); grid(headingAxes, 'on'); title(headingAxes, 'Heading'); ylim(headingAxes, [-90 90]); xlim(headingAxes, [t(1) t(end)]); sgtitle(fig, 'Targets Speed and Heading Estimates') for it = 1:numel(trackIDs) id = trackIDs(it); vx = targetStateEstimates(2, :, id); vy = targetStateEstimates(4, :, id); speedEstimate = sqrt(vx.^2 + vy.^2); headingEstimate = atan2d(vy, vx); plot(speedAxes, t, speedEstimate, 's--', 'LineWidth', 2, 'DisplayName', sprintf('Track %d estimate', trackIDs(it))); plot(headingAxes, t, headingEstimate, 's--', 'LineWidth', 2); end for it = 1:numel(trajectories) [~, ~, velocity] = lookupPose(trajectories{it}, t); trueSpeed = sqrt(velocity(:, 1).^2 + velocity(:, 2).^2); trueHeading = atan2d(velocity(:, 2), velocity(:, 1)); plot(speedAxes, t, trueSpeed, '-', 'LineWidth', 2, 'DisplayName', sprintf('Target %d truth', it)); plot(headingAxes, t, trueHeading, '-', 'LineWidth', 2); end end function R = rotationMatrix(angles) % Determine the "ZYX" axis-sequence rotation matrix from the set of % bearing, downtilt, and slant angles specifying the array orientation R = rotz(angles(1))*roty(angles(2))*rotx(angles(3)); end function ST = applyInitialPose(ST,trajectory) ST.Position = trajectory.Waypoints(1,:); ST.Orientation = [0 0 0]; ST.Velocity = trajectory.Velocities(1,:); end
