main.m-主函数文件（脚本文件用到的各个变量及响应赋值，和公式）
clear;
sample_csi_traceTmp =load('sample_csi_trace');//加载了名为 sample_csi_trace.mat 的文件，并将其存储在名为 sample_csi_traceTmp 的结构体变量(将.mat文件打开，送给结构体）
sample_csi_trace = sample_csi_traceTmp.sample_csi_trace;//文件中存储的变量名为 sample_csi_trace，所以通过 sample_csi_traceTmp.sample_csi_trace 来访问这个变量，并将其赋值给 sample_csi_trace 变量。（将文件中存储的东西通过结构体拿出来，赋值给一个变量，本质就是文件的内容无法直接送出，要经过结构体的打开）
fc = 5.63e9;//中心频率
M = 3;//rx天线个数
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子载波指数
N = length(SubCarrInd);//子载波数
fgap = 312.5e3;//WiFi中连续子载波之间以Hz为单位的频率差距
lambda = c/fc;//波长
T = 1;//发射机天线数
// MUSIC算法需要估计网格中的MUSIC频谱。paramRange捕获该网格的参数
//对于下面的示例，MUSIC频谱计算101ToF(飞行时间)值，间隔在- 25ns和25ns之间。MUSIC频谱是为101 AoA(到达角)值在-90度和90度之间计算的。
paramRange = struct;
paramRange.GridPts = [101 101 1];//格式格点数[ToF(飞行时间)格点数，到达角(AoA)格点数，1]
paramRange.delayRange = [-50 50]*1e-9;//考虑ToF网格的最低和最高值。
paramRange.angleRange = 90*[-1 1];//最低和值考虑的AoA网格。
do_second_iter = 0;//是否进行二次迭代。
paramRange.K = floor(M/2)+1;//平滑相关参数。
paramRange.L = floor(N/2);//平滑相关参数。
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重构成N*M的矩阵
[PhsSlope, PhsCons]=removePhsSlope(csi_plot,M,SubCarrInd,N);//调用removePhsSlope函数（消除斜率影响文件）
ToMult = exp(1i*(-PhsSlope*repmat(SubCarrInd(:),1,M) - PhsCons*ones(N,M)));//对CSI信道进行去斜率处理，ToMult表示对应子载波上的相位修正值
csi_plot = csi_plot.*ToMult;//信号CSI通过信道，ToF与AOA会对信号相位产生影响。
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 重构成一个列向量，并赋值给变量 sample_csi_trace_sanitized；即将3D矩阵N*M*T变成一个一维向量。
//估计到达角的MUSIC算法
//aoaEstimateMatrix是(nComps x 5)矩阵，其中nComps是环境中的路径数量。第一列是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))//调用backscatterEstimationMusic函数， 根据 backscatterEstimationMusic 函数的定义，返回值应该是 [Parameters, varargout]，其中 Parameters 是一个包含估计的参数的向量，varargout 是一个可选的输出变量。因此，aoaEstimateMatrix 可以是 Parameters 或者 varargout{1}，取决于调用函数时是否传递了额外的输出参数。如果调用时没有传递额外的参数，则 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;//使用非凸优化库Cvxgen来计算相位斜率和常数。如果不使用Cvxgen，则使用矩阵计算来计算相位斜率和常数。（本函数使用矩阵计算相位斜率和常数）
if ~useCvxgen //useCvxgen =0，所以下面的if语句执行
PhsRelFirstPac = unwrap(angle(vna_response));//angle函数为复数数组 z 的每个元素返回区[-π,π]中的相位角，unwrap函数：相位展开以确保相邻相位变化在[-pi, pi]的范围内，一列。（将vna_response矩阵的每个元素的相位展开，并将结果存储在PhsRelFirstPac矩阵中）
//下面的代码用于使用相位相对时
//第一个包
for antIdForPhs = 1:M //PhsRelFirstPac(1,antIdForPhs)表示第antIdForPhs个天线的第一个数据点的相位角度，而PhsRelFirstPac(1,1)则表示第一个天线的第一个数据点的相位角度。因此，PhsRelFirstPac(1,antIdForPhs) - PhsRelFirstPac(1,1)就是第antIdForPhs个天线的第一个数据点的相位角度与第一个天线的第一个数据点的相位角度之差。
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重复数组副本（将SubCarrlnd（：）矩阵复制成M*1矩阵模样），ones（N*M，1）生产一个N*M的全是1的矩阵。
x = A\PhsRelFirstPac(:);//x向量是A矩阵的伪逆与PhsRelFirstPac矩阵展开后的向量的乘积
PhsSlope = x(1);//x(1)是相位斜率
PhsCons = x(2);//x(2)是相位常数
end
backscatterEstimationMusic.m-MUSIC算法的回波功率谱密度估计（MUSIC用于从数据中估计信号频率和空间分量的算法，回波功率谱密度估计是分析信号的回波信号，即从接收到的信号中提取有用的信息）
function [Parameters, varargout] = backscatterEstimationMusic(csi_trace, M, N, c, fc, Ttot, fgap, SubCarrInd, d, paramRange, maxRapIters, useNoise, do_second_iter, varargin)
nComps = [];//初始化这个变量，将其设置为空值
corrThr = 0.4;//相关性阈值决定有多少组件后，我们决定切断MUSIC进程（决定在多少个成分后截断的相关阈值）
//以前使用的值为0.83。这个值应该是截止角的cosd与投影平面的平方
if ~isempty(varargin) //isempty判断矩阵是否为空，为空则返回逻辑1，否则0
parametersTrue = varargin{1};//将其赋值为varargin的第一个元素，{}表示cell数组，返回的是矩阵本身，即ones（2）
nComps = size(parametersTrue,1);//返回ones维度1的大小为2
if nComps==0, nComps = []; end;end
//音乐参数
//K =地板(M/2)+1;% K是选择用于平滑音乐的天线子集的数量。默认值为(M/2)+1
//L =地板(N/2);% L是为平滑音乐而选择的子载波子集的数目。默认值为floor(N/2)
//对于移动定位
if ~isfield(paramRange,'K')//判断K是否是结构体数组paramRange的一个字段的名称，是则返回1，不是则为0。
K = floor(M/2)+1;//y因为返回的是1，所以下面三行不执行
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是否是结构体数组paramRange的一个字段的名称，是则返回1，不是则为0。
delayRange = [-25 25]*1e-9;%[-25 70]//isfield因为返回的是1，所以下面三行不执行
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//是否进行两次迭代，因为否，所以直接执行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;
endelse
seconditerGridPts = [];//seconditerGridPts是一个空矩阵
end
if ~isfield(paramRange,'circularTx')//判断circularTx是否是结构体数组paramRange的一个字段的名称，是则返回1，不是则为0。这是0。
paramRange.circularTx = 0;
end
if paramRange.circularTx == 1//此行不执行，因为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);//调用gridVecBackscatter函数
EigDiffCutoff = 4;
if ~isfield(paramRange,'X')
X = formatCSI(csi_trace, N, M, Ttot, K, L, T);//调用formatCSI函数
else
X = paramRange.X;end
[Pn,Ps,Qn,Qs,EigenInfo]=GetQnBackscatter(X,EigDiffCutoff, nComps);//获取信号中反向散射信号的贡献,输出的结果包括噪声主子空间Pn，信号主子空间Ps，噪声子空间对应的投影矩阵Qn，信号子空间对应的投影矩阵 Qs，以及一些特征信息EigenInfo。
delayFromMusic =[];
angleFromMusic =[];
deltaFromMusic =[];
if ~useNoise
nIters = min(maxRapIters,size(Qs,2));
else
nIters = 1;end
maxCorr = zeros(nIters,1);
doBreak = 1;//是否应用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表示第一个输参数
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);
endend
alphaFromMusic = zeros(length(delayFromMusic),1);//创建一个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
endParameters = [delayFromMusic*1e9 angleFromMusic deltaFromMusic*1e3 abs(alphaFromMusic) angle(alphaFromMusic)];//将delayFromMusic *1e9 、angleFromMusic、deltaFromMusic*1e3 、abs(alphaFromMusic)、angle(alphaFromMusic)这五个变量沿着列方向拼接成一个矩阵，然后将这个矩阵作为Parameters变量的值
上文中包含了几个基础的函数：gridVecBackscatter用于计算回波信号在网格点处的分量值；formatCSI：对输入CSI数据进行格式化处理；GetQnBackscatter：获取Qn和Qs矩阵以进行MUSIC算法处理； gridSampleBackscatter：对估计结果进行后处理。
gridVecBackscatter-提供回波的模拟数据
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函数返回矩阵行列大小；ones返回矩阵1*3的全是1的矩阵；max将1与后面矩阵每一项对比，选最大的一个，返回1*3的矩阵；进行点除，返回最后的值[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);//将线性索引换成下标，本质就是将GridPts矩阵101*101*1矩阵排数，从上往下，再从左往右排数，再将从中选出1-numGridPoints，将其中第一位数送给delayIndices，第二位数送给angleIndices，第三位数送给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 = [];//将其设置为空矩阵
if ~generateAtot
load('aTotSaveSpotfi');else
aTot = zeros(N*M*T,numGridPoints);//组件一个90*（101*101）的全0数组。
for gridEntry = 1:numGridPoints
aTot(:,gridEntry) =gridSampleBackscatter(fc,T,deltaGridValue(gridEntry), M, u_sGridValue(gridEntry), c,SubCarrInd,fgap, delayGridValue(gridEntry));//调用gridSampleBackscatter函数
end
if generateAtot == 2
save('aTotSaveSpotfi','aTot')//将名为 aTot 的变量保存到名为 aTotSaveSpotfi 的 MATLAB 数据文件中；保存的文件格式通常为.mat。
end
end
gridSampleBackscatter-对接收到的信号进行波束形成，以便在目标检测和成像等应用中提取目标信息。（波束形成：通过对接收信号的加权和相位调节，从而实现指向性的信号辐射或接收）
function steeringVec = gridSampleBackscatter(fc, T, deltak, M, u_s, c, SubCarrInd, fgap, delay)
//fc:中心频率
//T:等位移的时间瞬间数
//delta:沿发射机tak方向的位移
aoaSteering = exp(-1i*2*pi*u_s*(-(M-1)/2:(M-1)/2)');//这行代码定义了一个名为 aoaSteering 的向量，该向量用于在入射角处生成一个指向性相应。
aoaSteeringInv = exp(-1i*2*pi*fc*(0:(T-1))'*deltak/c);//这行代码定义了一个名为 aoaSteeringInv 的向量，用于在信号采样的时间/频率上生成一个指向性响应。
delaySteering = exp(-1i*2*pi*(SubCarrInd(:))*fgap*delay);//这行代码定义了一个名为delaySteering 的向量，对应于在给定子载波下的延迟。
steeringVecDelayAoATmp = delaySteering*aoaSteering.';//这行代码定义了一个名为 steeringVecDelayAoATmp 的矩阵。
steeringVecDelayAoA = steeringVecDelayAoATmp(:);//这行代码将矩阵 steeringVecDelayAoATmp 压缩成一个列向量。
steeringVecTmp = steeringVecDelayAoA*aoaSteeringInv.';//这行代码定义了一个名为 steeringVecTmp 的矩阵，它是将 steeringVecDelayAoA 与 aoaSteeringInv的共轭转置相乘
steeringVec = steeringVecTmp(:);//将 steeringVecTmp 矩阵展平为列向量 steeringVec
formatCSI-对CSI数据进行格式化处理（将输入的CSI矩阵按时间，天线和子载波进行分块，并将这些分块转换成一个用于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个时间步对应的所有元素组成的二维矩阵，即ysenarr成为一个N*M的二维矩阵。
[N,M] = size(ysenarr);
D = [];
for m=1:M
D{m} = hankel(ysenarr(1:L,m),ysenarr(L:N,m));//用指定的列和行向量创建一个非对称 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));//设置一个全零矩阵
for start_idx = 1:K
tmp =cat(2,D{start_idx:start_idx+M-K});//串联矩阵，选择垂直串联矩阵。
De((start_idx-1)*L+(1:L),:) = tmp;//将temp放在De对应位置
end
X = [X; De];//将所有的 De 拼接成一个大的矩阵 X
X3D{iT} = De;//将 De 存储在 X3D 的第 iT 个元素中
end
//如果需要执行 3D MUSIC 算法，则将所有的 De 拼接成一个三维矩阵 X3D，并将 X3D 中的子矩阵按照时间顺序拼接成一个大的矩阵 X
X = zeros(K*L*T, (M - K + 1)*(N - L + 1)*(Ttot - T + 1));//将X设置为一个全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-实现了获取输入矩阵的后向散射(Qn)、信号空间投影矩阵(Ps)、噪声空间投影矩阵(Pn)、正向散射（Qs）、EigenInfo特征信息（包括U矩阵，信号空间的多路径数量，奇异值）。
function [Pn,Ps,Qn,Qs,EigenInfo] = GetQnBackscatter(X,EigDiffCutoff, nComps)
[Utmp,D] = eig(X*X');//计算矩阵XX'的特征分解，其中Utmp是XX'的特征向量矩阵，D是XX'的特征值矩阵
D = abs(D);//将矩阵D中的元素全部取绝对值,取绝对值是为了确保所有的特征值都是非负的。
[Dtmp,I] = sort(diag(D), 'descend');//将对角线元素从大到小排序，并且返回排序后的结果和排序后每个元素的原始下标。diag(D)提取矩阵D的对角线元素，形成一个列向量，descend表示降序。Dtmp保存排好序的对角线元素值，I保存排序前每个元素的下标。
D = diag(Dtmp);//diag 函数可以将一个向量转换成对角矩阵。
U = Utmp(:,I);//这行代码将矩阵 Utmp 的列按从大到小的顺序排列，并将结果存储在矩阵 U 中
//设置常数
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函数返回参数的对数值，mean返回向量的均值。nSegments和subarraySize分别是输入数据X的段数和子阵列大小。lambdaTot是信号子空间和噪声子空间的特征值向量。geo_mean(lambdaTot(k:end))/mean(lambdaTot(k:end))计算了特征值的几何平均值除以算术平均值。
end
[~, SignalEndIdxTmp] = min(MDL);//找到MDL向量最小值的索引，~表示忽略第一个输出参数，即最小值，SignalEndIdxTmp赋值为MDL的最小值的索引。
if useMDL
SignalEndIdx = max(SignalEndIdxTmp-1, 1);else
VecForSignalEnd = wkeep(diag(D),5,'l');//对角线元素向量 diag(D) 的长度限制为5，wkeep 函数的第二个参数表示保留的长度，第三个参数表示从左边开始保留元素。
diag(D(1:floor(length(D)/2)));//这行代码的作用是取矩阵D的前一半对角线元素，并将它们放在一个对角矩阵中。floor意思向下取整。
Criterion1 = diff(db(VecForSignalEnd))<=max(-EigDiffCutoff,min(diff(db(VecForSignalEnd))));//db是用于计算分贝 (dB) 的函数，其计算方式为:db（x）=20log10(x)；diff函数可以对向量进行差分操作，即对相邻元素之间的差值进行计算，输出一个长度比原向量小1的向量；表示满足一定条件的子空间信号成分的终止点的位置，
Criterion3 = (VecForSignalEnd(1:end-1)/VecForSignalEnd(1)>0.03);//元素的值是VecForSignalEnd中的每个元素除以VecForSignalEnd的第一个元素是否大于0.03，即对于VecForSignalEnd的每个元素，如果它除以VecForSignalEnd的第一个元素大于0.03，则相应的Criterion3元素为1，否则为0。
SignalEndIdx = find(Criterion1 & Criterion3,1,'last');//通过使用find函数，查找同时满足Criterion1和Criterion3条件的最后一个位置，并将该位置赋值给变量SignalEndIdx。通过同时满足这两个条件，可以找到信号的结束位置。
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
//各种输出的结果公式。
Qn = U(:,SignalEndIdx+1:end);//将变换矩阵U的第SignaEndldx+1到最后一列的所有列向量提取出来，组成矩阵Qn。
Pn = (Qn*Qn');//Qn与Qn的共轭转置或转置相乘
Qs = U(:,1:SignalEndIdx);//将变换矩阵U的第1到第SignalEndldx的所有列向量提取出来，组成矩阵Qs
Ps = (Qs*Qs');//Qs与Qs的共轭转置或转置相乘
EigenInfo = struct;//定义为结构体
EigenInfo.Umatrix = U;//结构体变量定义为U
EigenInfo.nMPs = SignalEndIdx;//结构体变量定义为SignalEndldx
EigenInfo.singularValuesdB = db(diag(D));//结构体变量定义为D的对角矩阵的分贝计算函数
RAPMusicGridMaxBackscatter-执行ICFD算法并计算所需的输出，如延迟，角度，多普勒偏移量，最大相关系数和MUSIC谱（ICFD是一种基于MUSIC算法的信号处理技术，包含两个步骤：第一步是对雷达回波信号进行波束形成，以提取信号的方向信息和延迟信息。第二步是对所得到的方向和延迟信息进行迭代估计和校正，以减少多径效应的影响。 ICFD技术的主要优点是，它可以提高目标信号的信噪比，并减少多径效应对目标信号的影响）
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：最大相关系数;music_spectrum_plot：MUSIC谱。
maxCorr = [];//将其设置为空矩阵
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);//创建一个方阵，大小为Qs的行数，列数与行数相同，like 参数指定了创建的矩阵的数据类型和储存方式，它使用了一个复数类型的2+1i 作为参考类型。（eye函数创建一个对角线上全是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共轭转置矩阵与PsA转置矩阵的乘积的 每一行总和的列向量。
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);
endif ~useNoise
[maxCorr,loopVarMax] = max(music_spectrum);//maxCorr是music_spectrum中的最大值, loopVarMax是maxCorr在music_spectrum中的索引位置，也就是maxCorr所在的频率下标
[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));//重构 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
本文章仅个人学习使用。
