FMCW雷达多运动目标检测

FMCW雷达多运动目标检测

去年的时候做了一个基于FMCW雷达多运动目标检测的matlab程序，对多运动目标检测，以图像的方式显示结果。效果还可以，能够同时实现高速多目标、大距离范围内的目标检测，程序的链接为（在mathworks网站的File Exchange里面）：
FMCW雷达多运动目标检测

下面是几个结果：


程序在matlab2018b下正常运行，需要Phased Array System Toolbox工具箱。
回波信号产生程序如下：

% This program generate the FMCW signal reflected from targets, mainly
% modified from the example of "phased.RangeEstimator System
% object" in matlab 2018b. It need the "Phased Array System Toolbox" to
% work.

clear all;
close all;
fs = 153.6e6;
c = physconst('LightSpeed');
M=256;
fc = c*M;
landa=c/fc;
prf=fs/1000;
pri=1/prf;
Np = 300;
RangeMax=c/2/prf;   %
VelMax=prf/2*landa/2;
DutyCycle=0.02;
Rresolution=c/2/(fs);

Numtgts = 5;
tgtpos = zeros(3,Numtgts);
tgtvel = zeros(3,Numtgts);
tgtpos(1,:) = [1 300 200 100 900];
tgtvel(1,:) = [110 0 -90 -80 60];
tgtrcs = db2pow(10)*[1 1 1 1 1];

tgtmotion = phased.Platform(tgtpos,tgtvel);
target = phased.RadarTarget('PropagationSpeed',c,'OperatingFrequency',fc, ...
    'MeanRCS',tgtrcs);
radarpos = [0;0;0];
radarvel = [0;0;0];
radarmotion = phased.Platform(radarpos,radarvel);

txantenna = phased.IsotropicAntennaElement;
rxantenna = clone(txantenna);

bw = fs/2;
waveform = phased.LinearFMWaveform('SampleRate',fs, ...
    'PRF',prf,'OutputFormat','Pulses','NumPulses',1,'SweepBandwidth',fs/2, ...
    'DurationSpecification','Duty cycle','DutyCycle',DutyCycle);
sig = waveform();
Nr = length(sig);
bwrms = bandwidth(waveform)/sqrt(12);
rngrms = c/bwrms;

peakpower = 10;
txgain = 36.0;
transmitter = phased.Transmitter( ...
    'PeakPower',peakpower, ...
    'Gain',txgain, ...
    'InUseOutputPort',true);
radiator = phased.Radiator( ...
    'Sensor',txantenna,...
    'PropagationSpeed',c,...
    'OperatingFrequency',fc);

channel = phased.FreeSpace( ...
    'SampleRate',fs, ...    
    'PropagationSpeed',c, ...
    'OperatingFrequency',fc, ...
    'TwoWayPropagation',true);

collector = phased.Collector( ...
    'Sensor',rxantenna, ...
    'PropagationSpeed',c, ...
    'OperatingFrequency',fc);
rxgain = 42.0;
noisefig = 1;
receiver = phased.ReceiverPreamp( ...
    'SampleRate',fs, ...
    'Gain',rxgain, ...
    'NoiseFigure',noisefig);

dt = pri;
cube = zeros(Nr,Np);
for n = 1:Np
    [sensorpos,sensorvel] = radarmotion(dt);
    [tgtpos,tgtvel] = tgtmotion(dt);
    [tgtrng,tgtang] = rangeangle(tgtpos,sensorpos);
    sig = waveform();
    [txsig,txstatus] = transmitter(sig);    %主要是发射增益控制
    txsig = radiator(txsig,tgtang);         %主要是经天线发射
    txsig = channel(txsig,sensorpos,tgtpos,sensorvel,tgtvel);    %信道
    tgtsig = target(txsig);     %经目标反射
    rxcol = collector(tgtsig,tgtang);
    rxsig = receiver(rxcol);    %主要是接收增益控制
    cube(:,n) = rxsig;
end
save fmcwmatlab2018data;
% imagesc([0:(Np-1)]*pri*1e6,[0:(Nr-1)]/fs*1e6,abs(cube))
% xlabel('Slow Time {\mu}s')
% ylabel('Fast Time {\mu}s')
% axis xy
% 
% ndop = 128;
% rangedopresp = phased.RangeDopplerResponse('SampleRate',fs, ...
%     'PropagationSpeed',c,'DopplerFFTLengthSource','Property', ...
%     'DopplerFFTLength',ndop,'DopplerOutput','Speed', ...
%     'OperatingFrequency',fc);
% matchingcoeff = getMatchedFilter(waveform);
% [rngdopresp,rnggrid,dopgrid] = rangedopresp(cube,matchingcoeff);
% figure;imagesc(dopgrid,rnggrid,10*log10(abs(rngdopresp)))
% xlabel('Closing Speed (m/s)')
% ylabel('Range (m)')
% axis xy
% 
% mfgain = matchingcoeff'*matchingcoeff;
% dopgain = Np;
% noisebw = fs;
% noisepower = noisepow(noisebw,receiver.NoiseFigure,receiver.ReferenceTemperature);
% noisepowerprc = mfgain*dopgain*noisepower;
% noise = noisepowerprc*ones(size(rngdopresp));
% 
% rangeestimator = phased.RangeEstimator('NumEstimatesSource','Auto', ...
%     'VarianceOutputPort',true,'NoisePowerSource','Input port', ...
%     'RMSResolution',rngrms);
% dopestimator = phased.DopplerEstimator('VarianceOutputPort',true, ...
%     'NoisePowerSource','Input port','NumPulses',Np);
% 
% detidx = NaN(2,Numtgts);
% tgtrng = rangeangle(tgtpos,radarpos);
% tgtspd = radialspeed(tgtpos,tgtvel,radarpos,radarvel);
% tgtdop = 2*speed2dop(tgtspd,c/fc);
% for m = 1:numel(tgtrng)
%     [~,iMin] = min(abs(rnggrid-tgtrng(m)));
%     detidx(1,m) = iMin;
%     [~,iMin] = min(abs(dopgrid-tgtspd(m)));
%     detidx(2,m) = iMin;
% end
% 
% ind = sub2ind(size(noise),detidx(1,:),detidx(2,:));
% 
% [rngest,rngvar] = rangeestimator(rngdopresp,rnggrid,detidx,noise(ind))
% [spdest,spdvar] = dopestimator(rngdopresp,dopgrid,detidx,noise(ind))

回波处理的函数如下：

clear all;
close all;
load fmcwmatlab2018data;
x=cube;
[nrn,nan]=size(x);
matchingcoeff=getMatchedFilter(waveform);
mlength=length(matchingcoeff);
lambda=c/fc;

fx=fftshift((linspace(-prf/2,prf/2-prf/nan,nan))*2*pi);
Vel=-fx/2*lambda/2/pi;
ts=(0/fs:1/fs:(nrn-1)/fs).';

% Eliminating the frequency shift in frequncy domain
x=fft(x,nrn,1);
x=fft(x,nan,2);
u=exp(1i*4*pi/lambda*(ts*Vel));  %
u=fft(u,nrn,1);
for k=1:nan
    x3=cconv(x(:,k),u(:,k),nrn);
    x(:,k)=x3/nrn;
end
x=ifft(x,nrn,1);
x=ifft(x,nan,2);

% Construct the signal generating matrix
z=zeros(nrn,nrn-mlength);
Rs=0:Rresolution:(nrn-1)*Rresolution;
L=(4*pi*Rs*2/lambda).^2;
L=L/max(abs(L));
L=1./L;
L(1)=L(2)*1;
L(end-mlength+1:end)=ones(mlength,1);
for k=1:nrn-mlength
    z(:,k)=circshift(sig,k-1)*L(k);
end
% the noise added to the signal
noiseadd=[wgn(nrn-mlength+1,mlength,20*log10(2.0324*1e-9),'complex');zeros(mlength-1,mlength)];
for k=1:mlength
    noiseadd(:,k)=circshift(noiseadd(:,k),k-1);
end
z=[z noiseadd];
% Solve the problem throught the Least Square method
x=pinv(z'*z)*z'*x;

% Because for close range and high speed targets, the side lobe in the velocity direction is very
% high, so the chebwin window is added to lower the side lobe. This step
% perhaps is  unnecessary for far range target such as farther than 5
% meters;
w=chebwin(nan);
for k=1:nrn
    x(k,:)=x(k,:).*w.';
end

x=fft(x,nan,2);
x=fftshift(x,2);

% Because for close range and high speed targets, the side lobe severely
% effect the detection, so for far range such as farther than 5 meters, these
% lines need not be zeroed.
x(:,1)=0*x(:,1);
x(:,2)=0*x(:,2);
x(:,end)=0*x(:,end);
x(:,end-1)=0*x(:,end-1);

V=fftshift(Vel);
[Vgrid,Rgrid]=meshgrid(V,Rs);
figure;mesh(Vgrid,Rgrid,abs(x));axis tight;
xlabel('Vel');ylabel('Range');
figure;imagesc(V,Rs,abs(x));grid on;axis tight;
xlabel('Vel');ylabel('Range');