雷达数据立方体表示方式

如上图所示，数据立方体的行对应一个脉冲的不同时间采样点（对应的于不同的距离），有多个采样点；数据立方体的列对应于接收机的通道（多个接收机，相控阵模式）；数据立方体的页对应于不同的脉冲。

二、建立数据立方体方法

1. 数据建模

1. Initialize a single radar target at a distance of 30 km from the radar with a velocity of 150 m/s
2. Create an eight-element array transmitter/receiver
3. Generate a linear FM pulse waveform with a sample rate of 1 MHz and PRF of 1 KHz
4. Describe the free space path loss of the propagation channel
5. Build the signal processing algorithms

The following code is used to set up the ULA and visualize the array:

% ULA Specs : 8 Element Uniform Linear Array with Cosine Antenna Element
antenna=phased.ULA;
antenna.NumElements = 8;   
cosineElement = phased.CosineAntennaElement;
antenna.Element = cosineElement;viewArray(antenna);
pattern(antenna,300e6,-180:180,0,...'Type','directivity',...'PropagationSpeed',3e8)

The code to set up the waveform is straightforward and concise:

% Waveform Specs
waveform=phased.LinearFMWaveform;
waveform.SampleRate = 1e6;                     % 1 MHz sample rate
waveform.PRF=1000;                             % Pulse Repetition Frequency of 1000
waveform.PulseWidth=1e-4;                      pulses per second  
nSamples = waveform.SampleRate/waveform.PRF;   % Define pulses on 1000 samples
y = step(waveform);
t = (0:nSamples-1)/waveform.SampleRate;

The code to set up the other system components used in the example is similarly concise, and each component contains all parameters needed to fully describe the section of the radar they represent.
The following loop generates the radar pulse, updates the target position, and propagates the pulse to the target and then back to the radar. All returned signals are collected using a receiver array, and the radar data cube is built up pulse by pulse in the loop. The step methods are used to generate or process the data for each component of the radar data cube.

%% Generate radar pulses 
for ii=1:nPulses 
wf=step(waveform); % Generate waveform 
[tgtPos, tgtVel] = step(PlatformModel,1/prf); % Update target position 
[~, tgtAng] = rangeangle(tgtPos, radarPos); % Calculate range/angle to target 
s0 = step(TX, wf); % Amplify signal 
s1 = step(txArray,s0, tgtAng); % Radiate the signal from the array 
s2 = step(ChannelModel, s1, radarPos, tgtPos, radarVel, tgtVel); % Propagate from radar to target and return 
s3 = step(TgtModel, s2); % Reflect signal from Target 
s4 = step(rxArray,s3,tgtAng); % Receive the signal at the array 
s5 = step(rxPreamp,s4); % Add rx noise 
datacube(:,:,ii) = s5(:,:); % Build data cube 1 pulse at a time 
end 

全流程

%% RadarCubePart1.m
% This program creates a radar datacube using functions from the Phased
% Array System Toolbox.  This program shows 2 examples:
%  1- A radar datacube with a single target, an 8 element uniform linear
%  array, 1000 samples per radar pulse, and 32 pulses
%  2- A radar datacube with 20 targets, a 121 element array mounted on the
%  surface of a sphere, 20k samples per pulse and 32 pulses
% In both cases the program is divided up into 2 sections. Section 1
% defines and initializes the objects needed to build the radar data cube,
% and Section 2 runs the simulation loop to generate the pulses, move the
% target(s), reflect the pulses off the target and back to the radar, 
% collect the received signals and build the datacube.% Copyright 2015, The MathWorks, Inc.% ULA Specs : 8 Element Uniform Linear Array with Cosine Antenna Element
antenna=phased.ULA;
antenna.NumElements = 8;   
cosineElement = phased.CosineAntennaElement;
antenna.Element = cosineElement;viewArray(antenna);
pattern(antenna,300e6,-180:180,0,...'Type','directivity',...'PropagationSpeed',3e8)%% -----*****************************-----------
pattern(antenna,300e6,[-180:180],[-90:90],'CoordinateSystem','polar', ...'Type','power',  'PropagationSpeed',3e8)
%% -----*****************************-----------
figure
pattern(antenna,300e6,[-180:180],[-90:90],'CoordinateSystem','polar')
%% ****************************************% Waveform Specs
waveform=phased.LinearFMWaveform;
waveform.SampleRate = 1e6;                     % 1 MHz sample rate
waveform.PRF=1000;                             % Pulse Repetition Frequency of 1000 pulses per second
waveform.PulseWidth=1e-4;
nSamples = waveform.SampleRate/waveform.PRF;   % Define pulses on 1000 samples
y = step(waveform);
t = (0:nSamples-1)/waveform.SampleRate;figure
plot(t,real(y));title('Radar Waveform'); xlabel('Time (sec)'); ylabel('Amplitude')
prf = waveform.PRF;% Transmitter Specs for amplification
TX=phased.Transmitter('Gain',20);% Target Specs
TgtModel=phased.RadarTarget;
tgtPos=[10e3*sqrt(3);10e3;0];              % Target at 20 km distance, 30 degree azimuth
tgtVel=[75*sqrt(3);75;0];                  % Radial velocity is 150 m/sec
% Platform Specs
PlatformModel=phased.Platform;
PlatformModel.InitialPosition = tgtPos;
PlatformModel.Velocity = tgtVel;% Channel Specs
ChannelModel = phased.FreeSpace;
ChannelModel.TwoWayPropagation=true;% Tx and Rx ArraytxArray = phased.Radiator(...'Sensor',antenna,...'OperatingFrequency',300e6);
rxArray = phased.Collector(...'Sensor',antenna,...'OperatingFrequency',300e6);
rxPreamp = phased.ReceiverPreamp('Gain',10,'NoiseFigure',5,...'SampleRate',1e6);% Variable definitions
radarPos = [0;0;0];
radarVel = [0;0;0];
nPulses = 32;
datacube = complex(zeros(nSamples,antenna.NumElements,nPulses));%% Generate radar pulses
for ii=1:nPulseswf=step(waveform);                                                 % Generate waveform[tgtPos, tgtVel] = step(PlatformModel,1/prf);                      % Update target position[~, tgtAng] = rangeangle(tgtPos, radarPos);                        % Calculate range/angle to targets0 = step(TX, wf);                                                 % Amplify signals1 = step(txArray,s0, tgtAng);                                     % Radiate the signal from the arrays2 = step(ChannelModel, s1, radarPos, tgtPos, radarVel, tgtVel);   % Propagate from radar to target and returns3 = step(TgtModel, s2);                                           % Reflect signal from Target s4 = step(rxArray,s3,tgtAng);                                      % Receive the signal at the arrays5 = step(rxPreamp,s4);                                            % Add rx noisedatacube(:,:,ii) = s5(:,:);                                        % Build data cube 1 pulse at a timeendt = (0:nPulses*nSamples-1)/waveform.SampleRate;
y = abs(datacube(:,1,:));
plot(t,y(:));title('Reflected Target Return (One Channel)'); xlabel('Time (sec)'); ylabel('Magnitude')%%

示例二：Simulate with 20 targets and 121 element array

clear all%% Simulate with 20 targets and 121 element array
nTgt = 20;          % 20 targets with random postions, velocities and sizes% Generate random positions, velocities and RCS values
tgtPos=[15e3+abs(randn(1,nTgt))*30e3; abs(randn(1,nTgt))*30e3;zeros(1,nTgt)];
tgtVel=[abs(randn(1,nTgt))*10e2;abs(randn(1,nTgt))*10e2;zeros(1,nTgt)];
RCS = abs(randn(20, 1))';TgtModel = phased.RadarTarget('MeanRCS',RCS);
Platform = phased.Platform('InitialPosition',tgtPos,'Velocity',tgtVel);% Conformal Array with 121 elements mounted on a spherical surface
az = -10:2:10;
el = -10:2:10;
[ele_az, ele_el] = meshgrid(az,el);
ele_az = ele_az(:).';
ele_el = ele_el(:).';
ele_normal = [ele_az; ele_el];
N_ele = size(ele_normal,2);
[ele_x, ele_y, ele_z] = sph2cart(degtorad(ele_az),degtorad(ele_el),ones(1,N_ele));
ele_pos = [ele_x;ele_y;ele_z]*10;% Display array geometry
antenna = phased.ConformalArray('ElementPosition',ele_pos,'ElementNormal',ele_normal);
viewArray(antenna)nElements = getNumElements(antenna);% Waveform Specs
waveform=phased.LinearFMWaveform;
waveform.SampleRate = 20e6;                       % 20 MHz sample rate
waveform.PRF=1000;                                % 20k samples per pulse
waveform.PulseWidth=1e-4;
prf = waveform.PRF;
nSamples = waveform.SampleRate/prf;% Tx and Rx Specs
TX=phased.Transmitter('Gain',20);txArray = phased.Radiator(...'Sensor',antenna);rxArray = phased.Collector(...'Sensor',antenna);rxPreamp = phased.ReceiverPreamp('Gain',10,'NoiseFigure',5,...'SampleRate',20e6);% Channel Specs
ChannelModel = phased.FreeSpace;
ChannelModel.TwoWayPropagation=true;
ChannelModel.SampleRate = waveform.SampleRate;% Variable definitions
radarPos = [0;0;0];
radarVel = [0;0;0];
nPulses = 32;
datacube = complex(zeros(nSamples,nElements,nPulses));
tgtAng = zeros(2,nTgt);
tgtVel = zeros(3,nTgt);
tgtRng = zeros(1,nTgt);
tic%% Generate radar pulses and assemble radar datacube
for ii=1:nPulseswf=step(waveform);                                                  % Generate waveform[tgtPos, tgtVel] = step(Platform,1/prf);                            % Update target positions[tgtRng, tgtAng] = rangeangle(tgtPos, radarPos);                    % Calculate range/angle to targetss0 = step(TX, wf);                                                  % Amplify signals1 = step(txArray,s0,tgtAng);                                       % Radiate the signal from Tx arrays2 = step(ChannelModel, s1, radarPos, tgtPos, radarVel, tgtVel);    % Propagate to target and returns3 = step(TgtModel, s2);                                            % Reflect signal from Target s4 = step(rxArray,s3,tgtAng);                                       % Receive signal at rx arrays5 = step(rxPreamp,s4);                                             % Add rx noisedatacube(:,:,ii) = s5(:,:);                                         % build data cube 1 pulse at a time
end
tocfigure 
t = (0:nPulses*nSamples-1)/waveform.SampleRate;
y = abs(datacube(:,1,:));
plot(t,y(:));title('Reflected Target Return (One Channel)'); xlabel('Time (sec)'); ylabel('Magnitude')figure
[t,r]=(cart2pol(tgtPos(1,:),tgtPos(2,:),tgtPos(3,:)));
polar(t,r,'x'); title('Target Positions')
%%
% Complete

三、利用数据立方体处理雷达数据

%% RadarCubePart2.m
% This program processes a radar datacube using functions from Phased Array
% System Toolbox.  The program first creates a radar datacube for a single
% target at 20 km distance, 30 degree azimuth and 150 m/sec velocity.  The
% program then performs beamforming, matched filtering and Doppler
% Processing to determine the target's distance and velocity.  % Copyright 2015, The MathWorks, Inc.% Target Specs
TgtModel=phased.RadarTarget;
tgtpos=[10e3*sqrt(3);10e3;0];  % Target at 20 km distance, 30 degree azimuth
tgtvel=[75*sqrt(3);75;0];      % Radial velocity is 150 m/sec% ULA Specs
antenna=phased.ULA;
antenna.NumElements = 8;
cosineElement = phased.CosineAntennaElement;
antenna.Element = cosineElement;% Waveform Specs
waveform=phased.LinearFMWaveform;
waveform.PRF=1000;
waveform.PulseWidth=1e-4;
prf = waveform.PRF;
nSamples = waveform.SampleRate/prf;% Transmitter Specs
TX=phased.Transmitter('Gain',20);% Platform Specs
PlatformModel=phased.Platform;
PlatformModel.InitialPosition = tgtpos;
PlatformModel.Velocity = tgtvel;% Channel Specs
ChannelModel = phased.FreeSpace;
ChannelModel.TwoWayPropagation=true;% Tx ans Rx SpecstxArray = phased.Radiator(...'Sensor',antenna,...'OperatingFrequency',300e6);
rxArray = phased.Collector(...'Sensor',antenna,...'OperatingFrequency',300e6);
rxPreamp = phased.ReceiverPreamp('Gain',10,'NoiseFigure',5);% Variable definitions
radarPos = [0;0;0];
radarVel = [0;0;0];nPulses = 32;
tgtAng = zeros(2,nPulses);
tgtAngcopy = zeros(2,nPulses);
datacube = complex(zeros(nSamples,antenna.NumElements,nPulses));%% Generate radar pulses
for ii=1:nPulseswf=step(waveform);                                                 % Generate waveform[tgtPos, tgtVel] = step(PlatformModel,1/prf);                      % Update target position[tgtRng, tgtAng] = rangeangle(tgtPos, radarPos);                   % Calculate range/angle to targettgtAngcopy(:,ii)=tgtAng;s0 = step(TX, wf);                                                 % Amplify signals1 = step(txArray,s0, tgtAng);                                     % Radiate the signal from the arrays2 = step(ChannelModel, s0, radarPos, tgtPos, radarVel, tgtVel);   % Propagate from radar to target and returns3 = step(TgtModel, s2);                                           % Reflect signal from Target s4 = step(rxArray,s3,tgtAng);                                      % Receive the signal at the arrays5 = step(rxPreamp,s4);                                            % Add rx noisedatacube(:,:,ii) = s5(:,:);                                        % Build data cube 1 pulse at a time
endfigure;t = (0:nPulses*nSamples-1)/waveform.SampleRate;
y = abs(datacube(:,1,:));
plot(t,y(:));title('Reflected Target Return (One Channel)'); xlabel('Time (sec)'); ylabel('Magnitude')%% Perform beamforming
% Beamformer Specs
beamformer=phased.PhaseShiftBeamformer;
beamformer.SensorArray=antenna;
beamformer.DirectionSource='Input port';
beamformer.WeightsOutputPort=true;
beamformer.WeightsNormalization='Preserve power';[bf0,w0]=step(beamformer,datacube(:,:,1),[0;0]);
[bf, w]=step(beamformer,datacube(:,:,1),[30;0]);figure;
subplot(2,2,1);
pattern(antenna,300e6,-180:180,0,...'PropagationSpeed',physconst('LightSpeed'),'Normalize',false,...'Type','powerdb','CoordinateSystem','rectangular');
subplot(2,2,2);plot(abs(bf0));
title('Sum Of Receive Elements'); xlabel('Time (msec)');
subplot(2,2,3);pattern(antenna,300e6,-180:180,0,'Weights',w,...'PropagationSpeed',physconst('LightSpeed'),'Normalize',false,...'Type','powerdb','CoordinateSystem','rectangular');
subplot(2,2,4);plot(abs(bf)); 
title('Beamformed Return (30 degree steering)'); xlabel('Time (msec)');%% Beamform for all 32 pulses
beamformed=complex(zeros(nSamples,nPulses));
for ii=1:nPulsesbeamformed(:,ii)=step(beamformer,squeeze(datacube(:,:,ii)),tgtAngcopy(:,ii));
end
figure; 
t = (0:nPulses*nSamples-1)/waveform.SampleRate;
y = abs(reshape(beamformed,nPulses*nSamples,1));
dc=abs(reshape(datacube(:,1,:),nPulses*nSamples,1));
subplot(2,1,1);plot(t,dc);title('Single Channel Target Return'); xlabel('Time (sec)'); ylabel('Magnitude')
subplot(2,1,2);plot(t,y);title('Beamformed Target Return'); xlabel('Time (sec)'); ylabel('Magnitude')%% Perform matched filtering
b = getMatchedFilter(waveform);matchedfilter = phased.MatchedFilter(...'Coefficients',b,...'SpectrumWindow','Hamming');matchFiltered = step(matchedfilter,beamformed);figure;
subplot(1,2,1);plot(real(wf));title('Radar Waveform'); xlabel('Time (sec)'); ylabel('Amplitude')
subplot(1,2,2);plot(real(b));title('Waveform Matched Filter');figure;
t = (0:nSamples-1)/waveform.SampleRate;
subplot(1,2,1);plot(t,abs(bf)); title('Beamformed Return (30 degree steering)')
subplot(1,2,2);plot(t,abs(matchFiltered(:,nPulses/2)));title('Pulse Compressed Return'); xlabel('Time (sec)'); ylabel('Magnitude')%% Find range bin where detected peak occurs
[m,ind] = max(abs(matchFiltered(:,nPulses/2)));  % ind is the range bin where the max amplitude occurs for the middle pulse
targetRange = time2range((ind-length(b)-1)/waveform.SampleRate, beamformer.PropagationSpeed)%% Perform doppler processing
dopplered = fftshift(fft(beamformed(ind,:).')); % Take the fft at the max amplitude range bin 
lambda=beamformer.PropagationSpeed/beamformer.OperatingFrequency;
h5 = figure;
f = (-prf/2:prf/nPulses:prf/2-prf/64);
v = f*lambda/2;
plot(v,abs(dopplered)); title('Doppler Processing'); xlabel('Target Speed (m/s)'); ylabel('Magnitude');
annotation(h5,'textbox',...[0.142657579062161 0.836671802773498 0.151671755725191 0.0647149460708782],...'String',strcat('Range = ',num2str(targetRange),'m'),...'FitBoxToText','on');%% Doppler Processing Using RangeDopplerResponse
% Use the RangeDopplerResopnse object to display target range and speed
figure;
rangeDoppler=phased.RangeDopplerResponse;
rangeDoppler.DopplerOutput='Speed';
plotResponse(rangeDoppler,beamformed,b)
ylabel('Range (km)');%% Doppler Processing Using AngleDopplerResponse
% First apply matched filter to each receive channel
matched = complex(zeros(nSamples,nPulses,antenna.NumElements));
for ii=1:antenna.NumElementsfor jj=1:nPulsesmatched(:,jj,ii) = filter(b,1,datacube(:,ii,jj));end
end% Now use the AngleDopplerResponse object to display Target Angle and 
% Dopppler frequency shift
dd = squeeze(matched(ind,:,:)).';
angleDopplerResp = phased.AngleDopplerResponse('SensorArray',antenna,...'OperatingFrequency',3e8, ...'PropagationSpeed',physconst('LightSpeed'),...'PRF',prf, 'ElevationAngle',0);
figure;
plotResponse(angleDopplerResp,dd);%%% Complete

四、相控阵运用示例

实现不同的绘图显示

c = 3e8;            %speed of light
range_max = 180;    %max detection range
tm = 6*(2*range_max/c); %sweep time
%tm is 7.2e-6 s
bw = 200e6;         %sweep bandwidth
sweep_slope = bw/tm;
v_max = 150*1000/3600;  %target max velocity
fc = 77e9;          %radar frequency
lambda = c/fc;      %radar wavelength
fs = 72e6;          %sampling rate
%sampling rate based on ADC datasheet
chirps = 64;        %frame size
samples = ceil(tm*fs);  %samples in one chirp%% target
R0 = 20; %range in meters
V = 40; %radial velocity, m/s%% 
t = 0; %time
mix = zeros(samples, chirps); %mixer output
for i=1:1:chirpstd = 2 * R0 / c; %round trip delayphi0 = 4*pi*fc*R0/c; %inital phaset = 0; % Resetfor j=1:1:samplesa = (-2*pi*fc*2*V*i*tm/c ...    %phase shift-2*pi*(2*V*(fc+i*bw)/c + sweep_slope*td)*t);   %frequencymix(j,i) = 0.5*cos(a);t = t + 1/fs;end
end%% Form the range-Doppler map (RDM)% RDM axes
rangeBinAxis = (0:samples-1).*c/(2*bw);
dopplerBinSize = (1/tm)/chirps;
velocityBinAxis = (-chirps/2:chirps/2-1).*dopplerBinSize*lambda/2;% 2D FFT to perform range and Doppler compression (i.e. form the RDM)
rdm = fftshift(fft2(mix), 2);% Plot the RDM for the valid ranges of interest - targets ahead of you
figure;
surf(velocityBinAxis, rangeBinAxis(1:ceil(samples/2)), 20*log10(abs(rdm(1:ceil(samples/2), :))));
% surf(velocityBinAxis, rangeBinAxis, 20*log10(abs(rdm)));  % See the entire spectrum
xlabel("Range (m)");
ylabel("Velocity (m/s)");
axis tight;
shading flat;
view(0, 90);
colorbar;% figure(1)
% rngdopresp = phased.RangeDopplerResponse('PropagationSpeed',c,...
%     'DopplerOutput','Speed','OperatingFrequency',fc,'SampleRate',fs,...
%     'RangeMethod','FFT','SweepSlope',sweep_slope,...
%     'RangeFFTLengthSource','Property','RangeFFTLength',2048,...
%     'DopplerFFTLengthSource','Property','DopplerFFTLength',256);
% 
% clf;
% plotResponse(rngdopresp,mix);
% axis([-v_max v_max 0 range_max])