main.m-主函數(shù)文件（腳本文件用到的各個變量及響應(yīng)賦值，和公式）

clear;

sample_csi_traceTmp =load('sample_csi_trace');//加載了名為 sample_csi_trace.mat 的文件，并將其存儲在名為 sample_csi_traceTmp 的結(jié)構(gòu)體變量(將.mat文件打開，送給結(jié)構(gòu)體）

sample_csi_trace = sample_csi_traceTmp.sample_csi_trace;//文件中存儲的變量名為 sample_csi_trace，所以通過 sample_csi_traceTmp.sample_csi_trace 來訪問這個變量，并將其賦值給 sample_csi_trace 變量。（將文件中存儲的東西通過結(jié)構(gòu)體拿出來，賦值給一個變量，本質(zhì)就是文件的內(nèi)容無法直接送出，要經(jīng)過結(jié)構(gòu)體的打開）

fc = 5.63e9;//中心頻率

M = 3;//rx天線個數(shù)

fs = 40e6;//通道頻寬

c = 3e8;//光速

d = 2.6e-2;//線性天線陣列中相鄰天線之間的距離

SubCarrInd =[-58,-54,-50,-46,-42,-38,-34,-30,-26,-22,-18,-14,-10,-6,-2,2,6,10,14,18,22,26,30,34,38,42,46,50,54,58];//CSI可用的WiFi子載波指數(shù)

N = length(SubCarrInd);//子載波數(shù)

fgap = 312.5e3;//WiFi中連續(xù)子載波之間以Hz為單位的頻率差距

lambda = c/fc;//波長

T = 1;//發(fā)射機天線數(shù)

//?MUSIC算法需要估計網(wǎng)格中的MUSIC頻譜。paramRange捕獲該網(wǎng)格的參數(shù)

//對于下面的示例，MUSIC頻譜計算101ToF(飛行時間)值，間隔在- 25ns和25ns之間。MUSIC頻譜是為101 AoA(到達(dá)角)值在-90度和90度之間計算的。

paramRange = struct;?

paramRange.GridPts = [101 101 1];//格式格點數(shù)[ToF(飛行時間)格點數(shù)，到達(dá)角(AoA)格點數(shù)，1]

paramRange.delayRange = [-50 50]*1e-9;//考慮ToF網(wǎng)格的最低和最高值。

paramRange.angleRange = 90*[-1 1];//最低和值考慮的AoA網(wǎng)格。

do_second_iter = 0;//是否進(jìn)行二次迭代。

paramRange.K = floor(M/2)+1;//平滑相關(guān)參數(shù)。

paramRange.L = floor(N/2);//平滑相關(guān)參數(shù)。

paramRange.T = 1;?

paramRange.deltaRange = [0 0];

maxRapIters = Inf;?

useNoise = 0;?

paramRange.generateAtot = 2;

//?ToF凈化代碼(SpotFi論文中算法1)

csi_plot = reshape(sample_csi_trace, N, M);//將sample_csi_trace重構(gòu)成N*M的矩陣

[PhsSlope, PhsCons]=removePhsSlope(csi_plot,M,SubCarrInd,N);//調(diào)用removePhsSlope函數(shù)（消除斜率影響文件）

ToMult = exp(1i*(-PhsSlope*repmat(SubCarrInd(:),1,M) - PhsCons*ones(N,M)));//對CSI信道進(jìn)行去斜率處理，ToMult表示對應(yīng)子載波上的相位修正值

csi_plot = csi_plot.*ToMult;//信號CSI通過信道，ToF與AOA會對信號相位產(chǎn)生影響。

relChannel_noSlope = reshape(csi_plot, N, M, T);//reshape用于改變矩陣/向量形狀；將csi_plot矩陣N*M變換成relChannel_noSlope矩陣N*M*T

sample_csi_trace_sanitized = relChannel_noSlope(:);//將矩陣relChannel_noSlope 重構(gòu)成一個列向量，并賦值給變量 sample_csi_trace_sanitized；即將3D矩陣N*M*T變成一個一維向量。

//估計到達(dá)角的MUSIC算法

//aoaEstimateMatrix是(nComps x 5)矩陣，其中nComps是環(huán)境中的路徑數(shù)量。第一列是ToF (ns)，第二列是AoA(度)，如SpotFi論文中定義的那樣

aoaEstimateMatrix = backscatterEstimationMusic(sample_csi_trace_sanitized, M, N, c, fc, T, fgap, SubCarrInd, d, paramRange, maxRapIters, useNoise, do_second_iter, ones(2))//調(diào)用backscatterEstimationMusic函數(shù)， 根據(jù) backscatterEstimationMusic 函數(shù)的定義，返回值應(yīng)該是 [Parameters, varargout]，其中 Parameters 是一個包含估計的參數(shù)的向量，varargout 是一個可選的輸出變量。因此，aoaEstimateMatrix 可以是 Parameters 或者 varargout{1}，取決于調(diào)用函數(shù)時是否傳遞了額外的輸出參數(shù)。如果調(diào)用時沒有傳遞額外的參數(shù)，則 aoaEstimateMatrix 是 Parameters，否則就是 varargout{1}。

tofEstimate=aoaEstimateMatrix(:,1);//從aoaEstimateMatrix中取出所有行的第一個元素，然后賦值給tofEstimate

aoaEstomate = aoaEstimateMatrix（:,2);//從 aoaEstimateMatrix 矩陣中取出第二列，即所有行的第二個元素，存儲到 aoaEstomate 向量中

removePhsSlope.m-消除相位斜率的影響文件

function[PhsSlope,PhsCons,vna_response_corrected]=removePhsSlope(vna_response,M,SubCarrInd,N)

useCvxgen = 0;//使用非凸優(yōu)化庫Cvxgen來計算相位斜率和常數(shù)。如果不使用Cvxgen，則使用矩陣計算來計算相位斜率和常數(shù)。（本函數(shù)使用矩陣計算相位斜率和常數(shù)）

if ~useCvxgen //useCvxgen =0，所以下面的if語句執(zhí)行

????PhsRelFirstPac = unwrap(angle(vna_response));//angle函數(shù)為復(fù)數(shù)數(shù)組?z?的每個元素返回區(qū)[-π,π]中的相位角，unwrap函數(shù)：相位展開以確保相鄰相位變化在[-pi, pi]的范圍內(nèi)，一列。（將vna_response矩陣的每個元素的相位展開，并將結(jié)果存儲在PhsRelFirstPac矩陣中）

//下面的代碼用于使用相位相對時

//第一個包

????for antIdForPhs = 1:M //PhsRelFirstPac(1,antIdForPhs)表示第antIdForPhs個天線的第一個數(shù)據(jù)點的相位角度，而PhsRelFirstPac(1,1)則表示第一個天線的第一個數(shù)據(jù)點的相位角度。因此，PhsRelFirstPac(1,antIdForPhs) - PhsRelFirstPac(1,1)就是第antIdForPhs個天線的第一個數(shù)據(jù)點的相位角度與第一個天線的第一個數(shù)據(jù)點的相位角度之差。

????????if PhsRelFirstPac(1,antIdForPhs) - PhsRelFirstPac(1,1) > pi?

????????????PhsRelFirstPac(:,antIdForPhs) = PhsRelFirstPac(:,antIdForPhs) - 2*pi;?

????????elseif PhsRelFirstPac(1,antIdForPhs) - PhsRelFirstPac(1,1) < -pi?

????????????PhsRelFirstPac(:,antIdForPhs) = PhsRelFirstPac(:,antIdForPhs) + 2*pi;?

????end

end?

A = [repmat(SubCarrInd(:), M, 1) ones(N*M, 1)]; //由所選的子載波編號和全1向量組成，repmat重復(fù)數(shù)組副本（將SubCarrlnd（：）矩陣復(fù)制成M*1矩陣模樣），ones（N*M，1）生產(chǎn)一個N*M的全是1的矩陣。

x = A\PhsRelFirstPac(:);//x向量是A矩陣的偽逆與PhsRelFirstPac矩陣展開后的向量的乘積

PhsSlope = x(1);//x(1)是相位斜率

PhsCons = x(2);//x(2)是相位常數(shù)

end

backscatterEstimationMusic.m-MUSIC算法的回波功率譜密度估計（MUSIC用于從數(shù)據(jù)中估計信號頻率和空間分量的算法，回波功率譜密度估計是分析信號的回波信號，即從接收到的信號中提取有用的信息）

function [Parameters, varargout] = backscatterEstimationMusic(csi_trace, M, N, c, fc, Ttot, fgap, SubCarrInd, d, paramRange, maxRapIters, useNoise, do_second_iter, varargin)

nComps = [];//初始化這個變量，將其設(shè)置為空值

corrThr = 0.4;//相關(guān)性閾值決定有多少組件后，我們決定切斷MUSIC進(jìn)程（決定在多少個成分后截斷的相關(guān)閾值）

//以前使用的值為0.83。這個值應(yīng)該是截止角的cosd與投影平面的平方

if ~isempty(varargin) //isempty判斷矩陣是否為空，為空則返回邏輯1，否則0

????parametersTrue = varargin{1};//將其賦值為varargin的第一個元素，{}表示cell數(shù)組，返回的是矩陣本身，即ones（2）

????nComps = size(parametersTrue,1);//返回ones維度1的大小為2

????if? ?nComps==0, nComps = []; end;

end

//音樂參數(shù)

//K =地板(M/2)+1;% K是選擇用于平滑音樂的天線子集的數(shù)量。默認(rèn)值為(M/2)+1

//L =地板(N/2);% L是為平滑音樂而選擇的子載波子集的數(shù)目。默認(rèn)值為floor(N/2)

//對于移動定位

if ~isfield(paramRange,'K')//判斷K是否是結(jié)構(gòu)體數(shù)組paramRange的一個字段的名稱，是則返回1，不是則為0。

????K = floor(M/2)+1;//y因為返回的是1，所以下面三行不執(zhí)行

????L = floor(N/2);?

????T = floor(Ttot/2)+1;

else?

????K = paramRange.K;?

????L = paramRange.L;?

????T = paramRange.T;?

end

//do_second_iter=0;//做兩次迭代

if ~isfield(paramRange,'delayRange')//判斷delayRange是否是結(jié)構(gòu)體數(shù)組paramRange的一個字段的名稱，是則返回1，不是則為0。

????delayRange = [-25 25]*1e-9;%[-25 70]//isfield因為返回的是1，所以下面三行不執(zhí)行

????deltaRange = [-c/2/fc c/2/fc];?

????angleRange = [-90 90];?

????GridPts = [100 100 100];??

else?

????delayRange = paramRange.delayRange;?

????deltaRange = paramRange.deltaRange;?

????angleRange = paramRange.angleRange;?

????GridPts = paramRange.GridPts;?

end?

MaxAngle = angleRange(2);//AOA相位角的最大值為90

MinAngle = angleRange(1);//AoA相位角的最小值為-90

if isempty(maxRapIters)//判斷maxRaplter是否為空，maxRaplters為Inf，則返回1.

?????maxRapIters = Inf;

end?

if do_second_iter//是否進(jìn)行兩次迭代，因為否，所以直接執(zhí)行else

????if~isfield(paramRange,'seconditerGridPts')?

????????GridPts = [70 70 35];

? ? ? ? MaxAngle = MaxAngle*(GridPts(2)-1)/(GridPts(2)+1);?

????????MinAngle = -MaxAngle;?

????????seconditerGridPts = [15 25 2]; %[15 25 5];?

????else?

????????seconditerGridPts = paramRange.seconditerGridPts;?

????end?

else?

????seconditerGridPts = [];//seconditerGridPts是一個空矩陣

?end?

if ~isfield(paramRange,'circularTx')//判斷circularTx是否是結(jié)構(gòu)體數(shù)組paramRange的一個字段的名稱，是則返回1，不是則為0。這是0。

????paramRange.circularTx = 0;?

end?

if paramRange.circularTx == 1//此行不執(zhí)行，因為paramRange.circularTx=0

????if ~isfield(paramRange, 'deltaRange')?

????????deltaRange = 0:359;?

????else?

????????deltaRange = paramRange.deltaRange;?

????end?

????dTx = paramRange.dTx;?

end?

[aTot,GridStart,GridSpacing, delayGridValue, u_sGridValue, deltaGridValue] = gridVecBackscatter(deltaRange, M, T, d, fc, c, do_second_iter, delayRange, SubCarrInd, N, fgap, GridPts, MaxAngle, MinAngle, paramRange.generateAtot);//調(diào)用gridVecBackscatter函數(shù)

EigDiffCutoff = 4;?

if ~isfield(paramRange,'X')

?????X = formatCSI(csi_trace, N, M, Ttot, K, L, T);//調(diào)用formatCSI函數(shù)

else?

????X = paramRange.X;?

end?

[Pn,Ps,Qn,Qs,EigenInfo]=GetQnBackscatter(X,EigDiffCutoff, nComps);//獲取信號中反向散射信號的貢獻(xiàn),輸出的結(jié)果包括噪聲主子空間Pn，信號主子空間Ps，噪聲子空間對應(yīng)的投影矩陣Qn，信號子空間對應(yīng)的投影矩陣 Qs，以及一些特征信息EigenInfo。

delayFromMusic =[];?

angleFromMusic =[];?

deltaFromMusic =[];

?if ~useNoise?

????nIters = min(maxRapIters,size(Qs,2));?

else?

????nIters = 1;?

end?

maxCorr = zeros(nIters,1);?

doBreak = 1;//是否應(yīng)用break語句

?if paramRange.circularTx == 0?

????for k=1:nIters?

????????[delayFromMusicTmp,angleFromMusicTmp, deltaFromMusicTmp, maxCorr(k),music_spectrum] =RAPMusicGridMaxBackscatter(aTot,GridStart,GridSpacing,GridPts,Qn,Qs,fc,fgap,d,K,L,delayFromMusic,angleFromMusic, deltaFromMusic,SubCarrInd, deltaGridValue, u_sGridValue, delayGridValue, T, c, do_second_iter, seconditerGridPts, useNoise);

? ? ? ? if k==1

????????????varargout{1} = music_spectrum;//varargout表示第一個輸參數(shù)

????????end?

????????if ~doBreak?

????????????delayFromMusic=[delayFromMusic; delayFromMusicTmp];?

????????????angleFromMusic=[angleFromMusic; angleFromMusicTmp];?

????????????deltaFromMusic=[deltaFromMusic; deltaFromMusicTmp];?

????????else?

????????????if k==1?

????????????????delayFromMusic = [delayFromMusic; delayFromMusicTmp];?

????????????????angleFromMusic = [angleFromMusic; angleFromMusicTmp];?

????????????????deltaFromMusic = [deltaFromMusic; deltaFromMusicTmp];?

????????????else?

????????????????if maxCorr(k)>corrThr*max(maxCorr)?

????????????????????delayFromMusic = [delayFromMusic; delayFromMusicTmp];?

????????????????????angleFromMusic = [angleFromMusic; angleFromMusicTmp];?

????????????????????deltaFromMusic = [deltaFromMusic; deltaFromMusicTmp];?

????????????????else?

????????????????????break?

????????????????end?

????????????end?

????????end?

????end?

????if ~doBreak?

????????allParameters = [delayFromMusic*1e9 angleFromMusic deltaFromMusic*1e3];?

????????varargout{1} = allParameters;?

????????cutoffEntry = find(maxCorr<corrThr*maxCorr(1),1)-1;?

????????if isempty(cutoffEntry)?

????????cutoffEntry = length(maxCorr);?

????????end?

????????delayFromMusic = delayFromMusic(1:cutoffEntry);?

????????angleFromMusic = angleFromMusic(1:cutoffEntry);?

????????deltaFromMusic = deltaFromMusic(1:cutoffEntry);

????end?

end?

alphaFromMusic = zeros(length(delayFromMusic),1);//創(chuàng)建一個length（delayFromMusic）*1的全零矩陣

if paramRange.circularTx == 1

????if ~isfield(paramRange, 'X')?

????????Ahat = [];?

????????for compNo = 1:length(delayFromMusic)?

????????u_s = (d*fc/c)*sind(angleFromMusic(compNo));?

????????Ahat(:,compNo) = gridSampleBackscatter(fc, Ttot, deltaFromMusic(compNo), M, u_s, c, SubCarrInd(1:N), fgap, delayFromMusic(compNo) );

?????????end?

????????alphaFromMusic = Ahat\csi_trace;?

????????residualEstimationError = norm(csi_trace - Ahat*alphaFromMusic)/norm(csi_trace);?

????????varargout{1}.residualEstimationError = residualEstimationError;?

????end?

end?

Parameters = [delayFromMusic*1e9 angleFromMusic deltaFromMusic*1e3 abs(alphaFromMusic) angle(alphaFromMusic)];//將delayFromMusic *1e9 、angleFromMusic、deltaFromMusic*1e3 、abs(alphaFromMusic)、angle(alphaFromMusic)這五個變量沿著列方向拼接成一個矩陣，然后將這個矩陣作為Parameters變量的值

上文中包含了幾個基礎(chǔ)的函數(shù)：gridVecBackscatter用于計算回波信號在網(wǎng)格點處的分量值；formatCSI：對輸入CSI數(shù)據(jù)進(jìn)行格式化處理；GetQnBackscatter：獲取Qn和Qs矩陣以進(jìn)行MUSIC算法處理； gridSampleBackscatter：對估計結(jié)果進(jìn)行后處理。

gridVecBackscatter-提供回波的模擬數(shù)據(jù)

function [aTot,GridStart,GridSpacing, delayGridValue, u_sGridValue, deltaGridValue] = gridVecBackscatter(deltaRange, M, T, d, fc, c, do_second_iter, delayRange, SubCarrInd, N, fgap, GridPts, MaxAngle, MinAngle, generateAtot)

GridStart = [delayRange(1) MinAngle deltaRange(1)];//輸入的值是[0,-90,0]

GridStop = [delayRange(2) MaxAngle deltaRange(2)];//輸入的值是[0,90,0]

GridSpacing = (GridStop-GridStart)./max(1, GridPts-ones(size(GridPts)));//size函數(shù)返回矩陣行列大小；ones返回矩陣1*3的全是1的矩陣；max將1與后面矩陣每一項對比，選最大的一個，返回1*3的矩陣；進(jìn)行點除，返回最后的值[0,1.8,0]。

angleStart = GridStart(2);//-90

delayStart = GridStart(1);//0

deltaStart = GridStart(3);//0

angleDiff = GridSpacing(2);//1.8

delayDiff = GridSpacing(1);//0

deltaDiff = GridSpacing(3);//0

numGridPoints = prod(GridPts);//prod返回向量元素乘積101*101*1

[delayIndices,angleIndices,deltaIndices] = ind2sub(GridPts,1:numGridPoints);//將線性索引換成下標(biāo)，本質(zhì)就是將GridPts矩陣101*101*1矩陣排數(shù)，從上往下，再從左往右排數(shù)，再將從中選出1-numGridPoints，將其中第一位數(shù)送給delayIndices，第二位數(shù)送給angleIndices，第三位數(shù)送給deltaIndices。

delayGridValue = delayStart + (delayIndices-1)*delayDiff;//0

angleGridValue = (angleStart + (angleIndices-1)*angleDiff)*pi/180;

u_sGridValue = (d*fc/c)*sin(angleGridValue);?

deltaGridValue = deltaStart + (deltaIndices-1)*deltaDiff;//0

aTot = [];//將其設(shè)置為空矩陣

if ~generateAtot

? ? load('aTotSaveSpotfi');

else

? ? aTot = zeros(N*M*T,numGridPoints);//組件一個90*（101*101）的全0數(shù)組。

? ? for gridEntry = 1:numGridPoints

? ? ? ? aTot(:,gridEntry) =gridSampleBackscatter(fc,T,deltaGridValue(gridEntry), M, u_sGridValue(gridEntry), c,SubCarrInd,fgap, delayGridValue(gridEntry));//調(diào)用gridSampleBackscatter函數(shù)

? ? end

? ? if generateAtot == 2

? ? ? ? save('aTotSaveSpotfi','aTot')//將名為 aTot 的變量保存到名為 aTotSaveSpotfi 的 MATLAB 數(shù)據(jù)文件中；保存的文件格式通常為.mat。

? ? end

end

gridSampleBackscatter-對接收到的信號進(jìn)行波束形成，以便在目標(biāo)檢測和成像等應(yīng)用中提取目標(biāo)信息。（波束形成：通過對接收信號的加權(quán)和相位調(diào)節(jié)，從而實現(xiàn)指向性的信號輻射或接收）

function steeringVec = gridSampleBackscatter(fc, T, deltak, M, u_s, c, SubCarrInd, fgap, delay)? ??

//fc:中心頻率

//T:等位移的時間瞬間數(shù)

//delta:沿發(fā)射機tak方向的位移

aoaSteering = exp(-1i*2*pi*u_s*(-(M-1)/2:(M-1)/2)');//這行代碼定義了一個名為 aoaSteering 的向量，該向量用于在入射角處生成一個指向性相應(yīng)。

aoaSteeringInv = exp(-1i*2*pi*fc*(0:(T-1))'*deltak/c);//這行代碼定義了一個名為 aoaSteeringInv 的向量，用于在信號采樣的時間/頻率上生成一個指向性響應(yīng)。

delaySteering = exp(-1i*2*pi*(SubCarrInd(:))*fgap*delay);//這行代碼定義了一個名為delaySteering 的向量，對應(yīng)于在給定子載波下的延遲。

steeringVecDelayAoATmp = delaySteering*aoaSteering.';//這行代碼定義了一個名為 steeringVecDelayAoATmp 的矩陣。

steeringVecDelayAoA = steeringVecDelayAoATmp(:);//這行代碼將矩陣 steeringVecDelayAoATmp 壓縮成一個列向量。

steeringVecTmp = steeringVecDelayAoA*aoaSteeringInv.';//這行代碼定義了一個名為 steeringVecTmp 的矩陣，它是將 steeringVecDelayAoA 與 aoaSteeringInv的共軛轉(zhuǎn)置相乘?

steeringVec = steeringVecTmp(:);//將 steeringVecTmp 矩陣展平為列向量 steeringVec

formatCSI-對CSI數(shù)據(jù)進(jìn)行格式化處理（將輸入的CSI矩陣按時間，天線和子載波進(jìn)行分塊，并將這些分塊轉(zhuǎn)換成一個用于3D MUSIC算法輸入的矩陣X）

function X = formatCSI(csi_trace, N, M, Ttot, K, L, T)

//X是我們輸入到MUSIC的矩陣

X = [];

X3D = [];//初始化兩個變量 X 和 X3D

csiTracePerPkt = reshape(csi_trace, N, M, Ttot);//并將輸入的 CSI 矩陣按照(N, M, Ttot)的形狀重新排列成 csiTracePerPkt

for iT = 1:Ttot

? ? ysenarr = csiTracePerPkt(:,:,iT);//選擇第iT個時間步對應(yīng)的所有元素組成的二維矩陣，即ysenarr成為一個N*M的二維矩陣。

? ? [N,M] = size(ysenarr);

? ? D = [];

? ? for m=1:M?

? ? ? ? D{m} = hankel(ysenarr(1:L,m),ysenarr(L:N,m));//用指定的列和行向量創(chuàng)建一個非對稱?Hankel矩陣,生成一個大小為 L x (N-L+1)的 Hankel 矩陣,對于 m=1:M 中的每個 m,D{m} 都存儲了一個 Hankel 矩陣。

? ? end

? ? De = zeros(K*L, (M - K +1)*(N - L+1));//設(shè)置一個全零矩陣

? ? for start_idx = 1:K

? ? ? ? tmp =cat(2,D{start_idx:start_idx+M-K});//串聯(lián)矩陣，選擇垂直串聯(lián)矩陣。

? ? ? ? De((start_idx-1)*L+(1:L),:) = tmp;//將temp放在De對應(yīng)位置

? ? end

? ? X = [X; De];//將所有的 De 拼接成一個大的矩陣 X

? ? X3D{iT} = De;//將 De 存儲在 X3D 的第 iT 個元素中

end

//如果需要執(zhí)行 3D MUSIC 算法，則將所有的 De 拼接成一個三維矩陣 X3D，并將 X3D 中的子矩陣按照時間順序拼接成一個大的矩陣 X

X = zeros(K*L*T, (M - K + 1)*(N - L + 1)*(Ttot - T + 1));//將X設(shè)置為一個全0矩陣

subarraySizeAntSubcarr = L*K;

for start_idx = 1:T

? ? tmp =? cat(2,X3D{start_idx:start_idx+Ttot-T});

? ? X((start_idx-1)*subarraySizeAntSubcarr + (1:subarraySizeAntSubcarr),:) = tmp;

end//最終輸出的矩陣 X 用于 MUSIC 算法的輸入。

GetQnBackscatter-實現(xiàn)了獲取輸入矩陣的后向散射(Qn)、信號空間投影矩陣(Ps)、噪聲空間投影矩陣(Pn)、正向散射（Qs）、EigenInfo特征信息（包括U矩陣，信號空間的多路徑數(shù)量，奇異值）。

function [Pn,Ps,Qn,Qs,EigenInfo] = GetQnBackscatter(X,EigDiffCutoff, nComps)

[Utmp,D] = eig(X*X');//計算矩陣XX'的特征分解，其中Utmp是XX'的特征向量矩陣，D是XX'的特征值矩陣

D = abs(D);//將矩陣D中的元素全部取絕對值,取絕對值是為了確保所有的特征值都是非負(fù)的。

[Dtmp,I] = sort(diag(D), 'descend');//將對角線元素從大到小排序，并且返回排序后的結(jié)果和排序后每個元素的原始下標(biāo)。diag(D)提取矩陣D的對角線元素，形成一個列向量，descend表示降序。Dtmp保存排好序的對角線元素值，I保存排序前每個元素的下標(biāo)。

D = diag(Dtmp);//diag 函數(shù)可以將一個向量轉(zhuǎn)換成對角矩陣。

U = Utmp(:,I);//這行代碼將矩陣 Utmp 的列按從大到小的順序排列，并將結(jié)果存儲在矩陣 U 中

//設(shè)置常數(shù)

minMP = 2;

useMDL = 0;

useDiffMaxVal = 0;

MDL = [];

lambdaTot = diag(D);//diag(D)提取矩陣D的對角線元素，形成一個列向量。

subarraySize = size(X,1);//返回X第一維度的大小。

nSegments = size(X,2);//返回X第二維度的大小。

maxMultipath = length(lambdaTot);//返回lanbdaTot的長度，即矩陣D對角線元素的長度。

for k = 1:maxMultipath

? ? MDL(k) =-nSegments*(subarraySize-(k-1))*log(geo_mean(lambdaTot(k:end))/mean(lambdaTot(k:end))) + 0.5*(k-1)*(2*subarraySize-k+1)*log(nSegments);//log函數(shù)返回參數(shù)的對數(shù)值，mean返回向量的均值。nSegments和subarraySize分別是輸入數(shù)據(jù)X的段數(shù)和子陣列大小。lambdaTot是信號子空間和噪聲子空間的特征值向量。geo_mean(lambdaTot(k:end))/mean(lambdaTot(k:end))計算了特征值的幾何平均值除以算術(shù)平均值。

end

[~, SignalEndIdxTmp] = min(MDL);//找到MDL向量最小值的索引，~表示忽略第一個輸出參數(shù)，即最小值，SignalEndIdxTmp賦值為MDL的最小值的索引。

if useMDL

? ? SignalEndIdx = max(SignalEndIdxTmp-1, 1);

else

? ? VecForSignalEnd = wkeep(diag(D),5,'l');//對角線元素向量 diag(D) 的長度限制為5，wkeep 函數(shù)的第二個參數(shù)表示保留的長度，第三個參數(shù)表示從左邊開始保留元素。

? ? diag(D(1:floor(length(D)/2)));//這行代碼的作用是取矩陣D的前一半對角線元素，并將它們放在一個對角矩陣中。floor意思向下取整。

? ? Criterion1 = diff(db(VecForSignalEnd))<=max(-EigDiffCutoff,min(diff(db(VecForSignalEnd))));//db是用于計算分貝 (dB) 的函數(shù)，其計算方式為:db（x）=20log10(x)；diff函數(shù)可以對向量進(jìn)行差分操作，即對相鄰元素之間的差值進(jìn)行計算，輸出一個長度比原向量小1的向量；表示滿足一定條件的子空間信號成分的終止點的位置，

? ? Criterion3 = (VecForSignalEnd(1:end-1)/VecForSignalEnd(1)>0.03);//元素的值是VecForSignalEnd中的每個元素除以VecForSignalEnd的第一個元素是否大于0.03，即對于VecForSignalEnd的每個元素，如果它除以VecForSignalEnd的第一個元素大于0.03，則相應(yīng)的Criterion3元素為1，否則為0。

? ? SignalEndIdx = find(Criterion1 & Criterion3,1,'last');//通過使用find函數(shù)，查找同時滿足Criterion1和Criterion3條件的最后一個位置，并將該位置賦值給變量SignalEndIdx。通過同時滿足這兩個條件，可以找到信號的結(jié)束位置。

? ? if isempty(SignalEndIdx)

? ? ? ? SignalEndIdx = find(Criterion3,1,'last');

? ? end

end

SignalEndIdx = max(SignalEndIdx, minMP);//minMP=2,找出SignalEndldx和minMP的最大值。

if ~isempty(nComps)//nComps=2,不是空

? ? SignalEndIdx = nComps;

end

//各種輸出的結(jié)果公式。

Qn = U(:,SignalEndIdx+1:end);//將變換矩陣U的第SignaEndldx+1到最后一列的所有列向量提取出來，組成矩陣Qn。

Pn = (Qn*Qn');//Qn與Qn的共軛轉(zhuǎn)置或轉(zhuǎn)置相乘

Qs = U(:,1:SignalEndIdx);//將變換矩陣U的第1到第SignalEndldx的所有列向量提取出來，組成矩陣Qs

Ps = (Qs*Qs');//Qs與Qs的共軛轉(zhuǎn)置或轉(zhuǎn)置相乘

EigenInfo = struct;//定義為結(jié)構(gòu)體

EigenInfo.Umatrix = U;//結(jié)構(gòu)體變量定義為U

EigenInfo.nMPs = SignalEndIdx;//結(jié)構(gòu)體變量定義為SignalEndldx

EigenInfo.singularValuesdB = db(diag(D));//結(jié)構(gòu)體變量定義為D的對角矩陣的分貝計算函數(shù)

RAPMusicGridMaxBackscatter-執(zhí)行ICFD算法并計算所需的輸出，如延遲，角度，多普勒偏移量，最大相關(guān)系數(shù)和MUSIC譜（ICFD是一種基于MUSIC算法的信號處理技術(shù)，包含兩個步驟：第一步是對雷達(dá)回波信號進(jìn)行波束形成，以提取信號的方向信息和延遲信息。第二步是對所得到的方向和延遲信息進(jìn)行迭代估計和校正，以減少多徑效應(yīng)的影響。?ICFD技術(shù)的主要優(yōu)點是，它可以提高目標(biāo)信號的信噪比，并減少多徑效應(yīng)對目標(biāo)信號的影響）

function [delayFromMusic, angleFromMusic, deltaFromMusic, maxCorr, music_spectrum_plot] = RAPMusicGridMaxBackscatter(aTot,GridStart,GridSpacing,GridPts,Qn,Qs,fc,fgap,d,K,L,delayFromMusicPrev,angleFromMusicPrev, deltaFromMusicPrev, SubCarrInd, deltaGridValue, u_sGridValue, delayGridValue, T, c, do_second_iter, seconditerGridPts, useNoise)//delayFromMusic：計算出的延遲向量;angleFromMusic：計算出的角度向量;deltaFromMusic：計算出的多普勒偏移量向量;maxCorr：最大相關(guān)系數(shù);music_spectrum_plot：MUSIC譜。

maxCorr = [];//將其設(shè)置為空矩陣

if ~useNoise

? ? if ~isempty(delayFromMusicPrev)

? ? ? ? Ahat = [];

? ? ? ? for compNo = 1:length(delayFromMusicPrev)

? ? ? ? ? ? u_s = (d*fc/c)*sind(angleFromMusicPrev(compNo));

? ? ? ? ? ? Ahat(:,compNo) = gridSampleBackscatter(fc, T, deltaFromMusicPrev(compNo), K, u_s, c, SubCarrInd(1:L), fgap, delayFromMusicPrev(compNo) );

? ? ? ? end

? ? ? ? PerpAhat = eye(size(Ahat,1)) - Ahat*pinv(Ahat'*Ahat)*Ahat';

? ? else

? ? ? ? PerpAhat = eye(size(Qs,1),'like',2+1i);//創(chuàng)建一個方陣，大小為Qs的行數(shù)，列數(shù)與行數(shù)相同，like 參數(shù)指定了創(chuàng)建的矩陣的數(shù)據(jù)類型和儲存方式，它使用了一個復(fù)數(shù)類型的2+1i 作為參考類型。（eye函數(shù)創(chuàng)建一個對角線上全是1，其他元素全是0的矩陣，然后乘以2+1i），最終得到的矩陣 PerpAhat 對角線上全是2+1i，其他元素全是0。

? ? end

? ? Qs = PerpAhat*Qs;

? ? Ps = (Qs*Qs');

? ? numGridPts = prod(GridPts);//返回向量元素乘積，即[101,101,1]乘積

? ? delayConsider = GridStart(1) + (0:GridPts(1)-1)*GridSpacing(1);

? ? u_sConsider = (d*fc/c)*sind( GridStart(2) + (0:GridPts(2)-1)*GridSpacing(2) );//sind返回元素的正弦

? ? deltaConsider = GridStart(3) + (0:GridPts(3)-1)*GridSpacing(3);

? ? delaySteeringMat = exp(-1i*2*pi*(SubCarrInd(1:L).')*fgap*delayConsider);

? ? aoaSteeringMat = exp(-1i*2*pi*((-(K-1)/2:(K-1)/2)')*u_sConsider);

? ? aoaSteeringInvMat = exp(-1i*2*pi*(fc/c)*((0:(T-1))')*deltaConsider);

? ? thetaTauMat = kron(aoaSteeringMat, delaySteeringMat);//返回矩陣aoaSteeringMat ?和 delaySteeringMat 的?Kronecker 張量積

? ? thetaTauDeltaMat = kron(aoaSteeringInvMat, thetaTauMat);

? ? PerpAhat_a = PerpAhat*thetaTauDeltaMat;

? ? PsA = Ps*PerpAhat_a;

? ? music_spec_num = sum((PerpAhat_a').*(PsA.'),2);//返回PerpAhat_a共軛轉(zhuǎn)置矩陣與PsA轉(zhuǎn)置矩陣的乘積的 每一行總和的列向量。

? ? music_spec_den = sum((PerpAhat_a').*(PerpAhat_a.'),2);

? ? music_spectrum = abs(music_spec_num./music_spec_den);//返回數(shù)組每個元素的絕對值

else

? ? delayConsider = GridStart(1) + (0:GridPts(1)-1)*GridSpacing(1);

? ? u_sConsider = (d*fc/c)*sind( GridStart(2) + (0:GridPts(2)-1)*GridSpacing(2) );

? ? deltaConsider = GridStart(3) + (0:GridPts(3)-1)*GridSpacing(3);

? ? delaySteeringMat = exp(-1i*2*pi*(SubCarrInd(1:L).')*fgap*delayConsider);

? ? aoaSteeringMat = exp(-1i*2*pi*((-(K-1)/2:(K-1)/2)')*u_sConsider);

? ? aoaSteeringInvMat = exp(-1i*2*pi*(fc/c)*((0:(T-1))')*deltaConsider);

? ? thetaTauMat = kron(aoaSteeringMat, delaySteeringMat);

? ? thetaTauDeltaMat = kron(aoaSteeringInvMat, thetaTauMat);

? ? PnA = (Qn*Qn')*thetaTauDeltaMat;

? ? thetaTauDeltaMatTrans = thetaTauDeltaMat';

? ? music_spectrum = sum(thetaTauDeltaMatTrans.*(PnA.'),2);

? ? music_spectrum = 1./abs(music_spectrum);

end

if ~useNoise

? ? [maxCorr,loopVarMax] = max(music_spectrum);//maxCorr是music_spectrum中的最大值, loopVarMax是maxCorr在music_spectrum中的索引位置，也就是maxCorr所在的頻率下標(biāo)

? ? [delay_idx,angle_idx, delta_idx] = ind2sub(GridPts,loopVarMax);//delay_idx 表示最大元素在第一維上的位置，angle_idx 表示最大元素在第二維上的位置，delta_idx 表示最大元素在第三維上的位置

? ? delayFromMusic = GridStart(1) + (delay_idx-1)*GridSpacing(1);

? ? angleFromMusic = GridStart(2) + (angle_idx-1)*GridSpacing(2);

? ? deltaFromMusic = GridStart(3) + (delta_idx-1)*GridSpacing(3);

else

? ? music_spectrum = reshape(music_spectrum,GridPts(1),GridPts(2),GridPts(3));

? ? BW = imregionalmax(music_spectrum);

? ? [delay_idx, angle_idx, delta_idx] = ind2sub(size(BW), find(BW));

? ? topPeakIndices = topk(music_spectrum(BW), min(size(Qs,2), length(find(BW(:)))));

? ? delayFromMusic = GridStart(1) + (delay_idx(topPeakIndices)-1)*GridSpacing(1);

? ? angleFromMusic = GridStart(2) + (angle_idx(topPeakIndices)-1)*GridSpacing(2);

? ? deltaFromMusic = GridStart(3) + (delta_idx(topPeakIndices)-1)*GridSpacing(3);

? ? maxCorr = max(music_spectrum);

end

music_spectrum = reshape(music_spectrum,GridPts(1),GridPts(2),GridPts(3));//重構(gòu) music_spectrum矩陣

music_spectrum_plot = music_spectrum;

if do_second_iter?

? ? delayFirstIter = delayFromMusic;

? ? angleFirstIter = angleFromMusic;

? ? deltaFirstIter = deltaFromMusic;

? ? for iComp = 1:length(delayFirstIter)

? ? ? ? GridPts = seconditerGridPts;

? ? ? ? delayRange = [delayFirstIter(iComp)-GridSpacing(1) delayFirstIter(iComp)+GridSpacing(1)];

? ? ? ? angleRange = [angleFirstIter(iComp)-GridSpacing(2) angleFirstIter(iComp)+GridSpacing(2)];

? ? ? ? deltaRange = [deltaFirstIter(iComp)-GridSpacing(3) deltaFirstIter(iComp)+GridSpacing(3)];

? ? ? ? [GridStart, GridSpacing, delayGridValue, u_sGridValue, deltaGridValue] = gridParamsBackscatter(GridPts, angleRange, deltaRange, d, fc, c, delayRange);

? ? ? ? if ~useNoise

? ? ? ? ? ? numGridPts = prod(GridPts);

? ? ? ? ? ? delayConsider = GridStart(1) + (0:GridPts(1)-1)*GridSpacing(1);

? ? ? ? ? ? u_sConsider = (d*fc/c)*sind( GridStart(2) + (0:GridPts(2)-1)*GridSpacing(2) );

? ? ? ? ? ? deltaConsider = GridStart(3) + (0:GridPts(3)-1)*GridSpacing(3);

? ? ? ? ? ? delaySteeringMat = exp(-1i*2*pi*(SubCarrInd(1:L).')*fgap*delayConsider);

? ? ? ? ? ? aoaSteeringMat = exp(-1i*2*pi*((-(K-1)/2:(K-1)/2)')*u_sConsider);

? ? ? ? ? ? aoaSteeringInvMat = exp(-1i*2*pi*(fc/c)*((0:(T-1))')*deltaConsider);

? ? ? ? ? ? thetaTauMat = kron(aoaSteeringMat, delaySteeringMat);

? ? ? ? ? ? thetaTauDeltaMat = kron(aoaSteeringInvMat, thetaTauMat);

? ? ? ? ? ? PerpAhat_a = PerpAhat*thetaTauDeltaMat;?

? ? ? ? ? ? PsA = Ps*PerpAhat_a;

? ? ? ? ? ? music_spec_num = sum((PerpAhat_a').*(PsA.'),2);? ? ? ? ? ??

? ? ? ? ? ? music_spec_den = sum((PerpAhat_a').*(PerpAhat_a.'),2);

? ? ? ? ? ? music_spectrum = abs(music_spec_num./music_spec_den);

? ? ? ?else

? ? ? ? ? ? delayConsider = GridStart(1) + (0:GridPts(1)-1)*GridSpacing(1);

? ? ? ? ? ? u_sConsider = (d*fc/c)*sind( GridStart(2) + (0:GridPts(2)-1)*GridSpacing(2) );

? ? ? ? ? ? deltaConsider = GridStart(3) + (0:GridPts(3)-1)*GridSpacing(3);

? ? ? ? ? ? delaySteeringMat = exp(-1i*2*pi*(SubCarrInd(1:L).')*fgap*delayConsider);

? ? ? ? ? ? aoaSteeringMat = exp(-1i*2*pi*((-(K-1)/2:(K-1)/2)')*u_sConsider);

? ? ? ? ? ? aoaSteeringInvMat = exp(-1i*2*pi*(fc/c)*((0:(T-1))')*deltaConsider);

? ? ? ? ? ? thetaTauMat = kron(aoaSteeringMat, delaySteeringMat);

? ? ? ? ? ? thetaTauDeltaMat = kron(aoaSteeringInvMat, thetaTauMat);

? ? ? ? ? ? PnA = (Qn*Qn')*thetaTauDeltaMat;

? ? ? ? ? ? thetaTauDeltaMatTrans = thetaTauDeltaMat';

? ? ? ? ? ? music_spectrum = sum(thetaTauDeltaMatTrans.*(PnA.'),2);

? ? ? ? ? ? music_spectrum = 1./abs(music_spectrum);

? ? ? ? end

? ? ? ? if ~useNoise

? ? ? ? ? ? [maxCorrTmp,loopVarMaxTmp] = max(music_spectrum);

? ? ? ? ? ? maxCorr = maxCorrTmp(1);

? ? ? ? ? ? loopVarMax = loopVarMaxTmp(1);

? ? ? ? ? ? [delay_idx,angle_idx, delta_idx] = ind2sub(GridPts,loopVarMax);

? ? ? ? ? ? delayFromMusic = GridStart(1) + (delay_idx-1)*GridSpacing(1);

? ? ? ? ? ? angleFromMusic = GridStart(2) + (angle_idx-1)*GridSpacing(2);

? ? ? ? ? ? deltaFromMusic = GridStart(3) + (delta_idx-1)*GridSpacing(3);

? ? ? ? else

? ? ? ? ? ? maxCorr = max(music_spectrum);

? ? ? ? ? ? music_spectrum = reshape(music_spectrum,GridPts(1),GridPts(2),GridPts(3));

? ? ? ? ? ? BW = imregionalmax(music_spectrum);

? ? ? ? ? ? [delay_idx, angle_idx, delta_idx] = ind2sub(size(BW), find(BW));

? ? ? ? ? ? topPeakIndices = topk(music_spectrum(BW), 1);

? ? ? ? ? ? delayFromMusic(iComp) = GridStart(1) + (delay_idx(topPeakIndices)-1)*GridSpacing(1);

? ? ? ? ? ? angleFromMusic(iComp) = GridStart(2) + (angle_idx(topPeakIndices)-1)*GridSpacing(2);

? ? ? ? ? ? deltaFromMusic(iComp) = GridStart(3) + (delta_idx(topPeakIndices)-1)*GridSpacing(3);

? ? ? ? end

? ? end

end

本文章僅個人學(xué)習(xí)使用。