logoChat01.ai

下面的框架中是一个OFDM ISAC的框架，主要关注delay和doppler的感知。下面写的是主函...

Criado em: 9 de janeiro de 2026

Respondido usando GPT-5.2 Thinking por Chat01

Pergunta

下面的框架中是一个OFDM ISAC的框架，主要关注delay和doppler的感知。下面写的是主函数。现在我们想改成afdm，并在ofdmmodulate和demodulate里面把fft换成daft，同时在接收端信道估计的时候也做了如下处理，利用2d-dft来估计我们的位置和速度。先看一下下面的代码并学习一下。
主函数：clc;
clear;
tic
addpath(genpath(pwd)) % include all the subfolders under the current folder into the calling function directory
snrdB = 10; % SNR configuration
simulationNum = 50;
for SNR = snrdB
rng('default');
[sensing] = sensingInit(); % sensing parameter initialization

text
[TxFreqSignal,sensing]=RsMapping(sensing); % sensing reference signal frequency domain mapping %% whether comm

rangeErrSet=zeros(sensing.targetNum,simulationNum);
veloErrSet=zeros(sensing.targetNum,simulationNum);
angleErrSet=zeros(sensing.targetNum,simulationNum);

for num=1:simulationNum
    [sensing,channel]=channelParaConfig(sensing,num); % channel parameter configuration %%target set

    [sensing,channel] = channelInit(sensing,channel); % initialization of channel parameters %% target add H

    [TxFreqSignalprecode] = TxAnalogPrecoding(sensing, TxFreqSignal, channel); % multi anteana care

    [TxTimeSignal,sensing]=ofdmModulate(TxFreqSignalprecode,sensing,num); % ofdm modulation

    channelOut=channelPass(TxTimeSignal,TxFreqSignalprecode,channel,sensing); % pass channel %%(time/frec)

    channelOut=awgn(channelOut,SNR); % add Gaussian noise

    [RxFreqSignal]=ofdmDemodulate(channelOut,sensing); % ofdm demodulation

    Hest=channelEst(TxFreqSignal,RxFreqSignal,sensing); % frequency domain channel estimation   %% (correlation)

    Hest=clutterCancel(Hest,sensing);  % Clutter cancellation
    
    [range,velocity,angle]=seningProcess(Hest,sensing); % sensing algorithm processing

    [rangeErr,veloErr,angleErr]=errStatistics(range,angle,velocity,sensing,channel);  % result statistics
    rangeErrSet(:,num)=rangeErr;
    veloErrSet(:,num)=veloErr;
    angleErrSet(:,num)=angleErr;

    % [Xposi,Yposi]=positionCalcu(range,angle,sensing,channel);
end
plotfig(sensing,rangeErrSet,veloErrSet,angleErrSet,channel,SNR); % draw the result figure

end
toc

function [Range,Velocity,Angle]=dft_2d(Hest,sensing)
targetNum=sensing.targetNum;
fc=sensing.fc;
Deltaf=sensing.SCS;
slotT=sensing.frameLen*1e-3/sensing.sensingSymbolNum;
c=3e8;
UEdV=sensing.UEdV;
UEdH=sensing.UEdH;
freqSpace=sensing.freqSpace;
pseudoPeakFlag=sensing.pseudoPeakFlag;
Range=zeros(targetNum,1);
Velocity=zeros(targetNum,1);
Angle=zeros(targetNum,1);
ang=0;

if sensing.overSamplFlag==1
lenf=sensing.overSamplNumf;
lent=sensing.overSamplNumt;
else
[lenf,lent]=size(Hest(:,:,1));
end
rxAntNum=sensing.rxAntNum;
initH=zeros(lenf,lent);
dataSet=cell(1,rxAntNum);
for antId=1:rxAntNum
Htemp=Hest(:,:,antId);
Hdft=fft(ifft(Htemp,lenf).',lent).';
dataSet{antId}=fftshift(Hdft,2);
initH=initH+fftshift(abs(Hdft),2);
end
residualH=initH;

for Idx=1:targetNum
loop=1;
flag=0;
while loop
if Idx>1 || flag==1
[Nc,Nt,Np]=size(Hest);

text
        for antId=1:rxAntNum
            for delayIdx=delayTapSet-lent/2-1
                for dolpIdx=dolpTapSet-lent/2-1
                    distance=delayIdx*deltaR;
                    phiA = exp(-1j*2*pi*[0:Np-1]'*0.5*cos(ang*pi/180));
                    phiR = exp(-1j*2*pi*Deltaf*freqSpace*distance*[0:Nc-1]/c);
                    phiV=exp(1j*2*pi*dolpIdx*2*deltafd*slotT*[0:Nt-1]);
                    phi=kron(phiR,phiV);
                    phiphi=kron(phiA',phi);
                    phi=phiphi;
                    reSignal=abs(fftshift(fft(ifft(phi,lenf).',lent).',2));
                    reSignal=reSignal/max(reSignal(:));
                    factor=residualH(delayIdx+1,dolpIdx+1+lent/2+1);
                    DDRxH = factor*reSignal;
                    residualH=residualH-abs(DDRxH);
                    residualH(residualH<0)=0;
                end
            end
        end
    end
end

% [Nc,Nt]=size(residualH);
% Vset= c*[-Nt/2:Nt/2-1]./(fc*slotT*Nt);
% Rset=[0:Nc-1]./(Deltaf*Nc)*c;
deltaV=c/(fc*slotT*lent*2);
deltafd=1/(lenf*slotT*2);
deltaR=c/(freqSpace*Deltaf*lenf);

% residualH=residualH(1:floor(lenf/2),:);
residualH=residualH(1:floor(lenf/1),:);
[~,index]=max(residualH(:));
[maxRow,maxCol]=ind2sub(size(residualH),index);

surf(residualH);
figure
plot(residualH(1:maxCol));

distance=(maxRow-1)*deltaR;
rangeReal=sensing.targetRange;
realRow = rangeReal/deltaR+1;

speed=(maxCol-lent/2-1)*deltaV;
if pseudoPeakFlag==1
    dolpTapSet=[maxCol-1,maxCol,maxCol+1];
    dolpTapSet(dolpTapSet>lent)=[];
    dolpTapSet(dolpTapSet<=0)=[];
    delayTapSet=[maxRow-1,maxRow,maxRow+1];
    delayTapSet(delayTapSet<=0)=[];
else
    dolpTapSet=maxCol;
    delayTapSet=maxRow;
end

% INN=120;
% mesh(deltafd*[-lent/2:lent/2-1],deltaR*[0:Nc-1],residualH(1:INN,:));
% xlim([-1000,1000])
% xlabel('Doppler[Hz]');
% ylabel('Range[m]');
% zlabel('Amplitude');

% angle estimate with MUSIC algorithm
if rxAntNum>1 && abs(speed)>=min(deltaV,5)
    Xset=zeros(1,rxAntNum);
    for rxAntId=1:rxAntNum
        Xset(1,rxAntId)=dataSet{rxAntId}(maxRow,maxCol);
    end
    RxT=Xset; % 补充Xset到RxT的赋值（原代码隐含此步骤）
    J=flip(eye(rxAntNum));
    clear Xset
    Rv=(RxT*conj(RxT)')/2;
    [U,Value]=eig(Rv);
    % Value=diag(Value);
    sigNum=targetNum;
    Un=U(:,1:end-sigNum);
    Pn=Un*Un';

    n=0;
    if sensing.verticalFlag==1
        Aset=sensing.AsetV;
        Psang=zeros(1,length(Aset));
        for A=Aset
            n=n+1;
            phiA = exp(1j*2*pi*[0:rxAntNum-1]*UEdV*cos(A*pi/180)); % ZOA estimation
            Psang(n)=1/(abs(phiA'*Pn*phiA));
        end
    else
        Aset=sensing.AsetH;
        Psang=zeros(1,length(Aset));
        for A=Aset
            n=n+1;
            phiA = exp(1j*2*pi*[0:rxAntNum-1]*UEdH*sin(A*pi/180)); % AOA estimation
            Psang(n)=1/(abs(phiA'*Pn*phiA));
        end
    end

    [peakPower,peakInd]=findpeaks(Psang);
    [~,index]=sort(peakPower,'descend');
    ang=Aset(peakInd(index(1:1)));
end

% if abs(speed)<min(deltaV,5)
    loop=loop+1;
    flag=1;
% else
    Range(Idx)=distance;
    Velocity(Idx)=speed;
    Angle(Idx)=ang;
    loop=0;
    flag=0;
end

if loop>10
    break;
end

end
end

function [TxTimeSignal,sensing]=ofdmModulate(TxFreqSignal,sensing,frameId)
Nsybmol=sensing.Nsymbol;
txAntNum=sensing.txAntNum;
fftSize=sensing.fftSize;
if sensing.blockCP==1
CplenSet=zeros(1,Nsybmol);
CplenSet(1:14:end)=fftSize;
sensing.CplenSet=CplenSet;
else
CplenSet=sensing.CplenSet;
end

removeCpInd=zeros(1,Nsybmol); % CP position index
TxTimeSignal=zeros(txAntNum,Nsybmol*fftSize+sum(CplenSet(1:Nsybmol)));

for antId=1:txAntNum
TxTimeSignalTemp=[];
% changed:choose AFD
% sensing.useAFDM
if sensing.useAFDM
TxSignalT=idaft(TxFreqSignal(:,:,antId),sensing.afdm.c1,sensing.afdm.c2); % idaft
else
TxSignalT=ifft(TxFreqSignal(:,:,antId),fftSize)*sqrt(fftSize); % ifft
end
%%%
for symInd=1:Nsybmol
removeCpInd(symInd)=length(TxTimeSignalTemp);
Cpdata=TxSignalT(end-CplenSet(symInd)+1:end,symInd);
TxSignalCp=[Cpdata.',TxSignalT(:,symInd).']; % add CP
TxTimeSignalTemp=cat(2,TxTimeSignalTemp,TxSignalCp);
end
TxTimeSignal(antId,:)=TxTimeSignalTemp;
end
removeCpInd(removeCpInd+1)=length(TxTimeSignalTemp);
sensing.removeCpInd=removeCpInd;
sensing.frameId=frameId;
end

function [RxFreqSignal]=ofdmDemodulate(channelOut,sensing)
if sensing.cfr==0
Nsybmol=sensing.Nsymbol;
fftSize=sensing.fftSize;
CplenSet=sensing.CplenSet;
removeCpInd=sensing.removeCpInd;
rxAntNum=sensing.rxAntNum;

text
% remove CP and FFT
RxSignal=zeros(fftSize,Nsybmol,rxAntNum);
RxFreqSignal=zeros(fftSize,Nsybmol,rxAntNum);

for antId=1:rxAntNum
    for symInd=1:Nsybmol
        RxSignal(:,symInd,antId)=channelOut(antId, removeCpInd(symInd)+CplenSet(symInd)+1 : removeCpInd(symInd)+CplenSet(symInd)+fftSize);
        %%% AFD added
        if sensing.useAFDM
            RxFreqSignal(:,symInd,antId)=daft(RxSignal(:,symInd,antId),sensing.afdm.c1,sensing.afdm.c2); % daft
        else
            RxFreqSignal(:,symInd,antId)=fft(RxSignal(:,symInd,antId))/sqrt(fftSize); % fft
        end
    end
end

else
RxFreqSignal=channelOut; % if pass frequency domain channel,skip this step
end
end
function Hest=channelEst(TxFreqSignal,RxFreqSignal,sensing)
timeSpace = sensing.timeSpace;
freqSpace = sensing.freqSpace;
Nzc=sensing.Nzc;
Nc=sensing.Nc;
fpattern=sensing.fpattern;
tpattern=sensing.tpattern;
sensingSymbolNum=sensing.sensingSymbolNum;
if sensing.fastchannel==1
tpattern=[1:sensingSymbolNum];
end

rxAntNum=sensing.rxAntNum;

%% AFD
if sensing.useAFDM == 1
txTime = idaft(TxFreqSignal,sensing.afdm.c1,sensing.afdm.c2);
rxTime = zeros(size(RxFreqSignal));
for rx=1:rxAntNum
rxTime(:,:,rx) = idaft(RxFreqSignal(:,:,rx),sensing.afdm.c1,sensing.afdm.c2);
end
TxFreqSignal = fft(txTime,Nc,1)/sqrt(Nc);
RxFreqSignal = fft(rxTime,Nc,1)/sqrt(Nc);
end

TxFreqSignal=circshift(TxFreqSignal,Nc/2,1);
TxSignalTF=TxFreqSignal(fpattern,tpattern);
Hest=zeros(length(sensingSymbolNum),sensingSymbolNum,rxAntNum);

RxFreqSignal=circshift(RxFreqSignal,Nc/2,1); % Frequency domain data passing through the time domain channel

for RxIdx=1:rxAntNum
RxSignalTF=RxFreqSignal(fpattern,tpattern,RxIdx);
Hest(:,:,RxIdx)=RxSignalTF ./ TxSignalTF;
end

if sensing.otfsSensingFlag==1
Hest=RxFreqSignal(1:Nc,tpattern);
end
end

Compartilhar esta Q&A