how can I create a simple clutter model?

what is the easy way to model clutter? without any toolbox. is there any mathematical expression for the clutter model. I want to see these radar signals effected from these clutters.

thanks a lot, have a good day.

Answer 1

[William Rose]

@Kerem Yüksel,

[edit: fix spelling errors]

[Edit 2: My initial implementation had an exponential distribution of clutter radii. I had intended for the areas to be exponentially distributed. That is fixed now.]

Yes there are models of radar clutter that do not require a toolbox.

Here is a good place to start:

Radar Clutter Modeling and Analysis (book chapter, 2014)

https://doi.org/10.1016/B978-0-12-396500-4.00011-9

I do not have immediate access to the full text above. Maybe you can get it through your university or employer.

The clutter model depends on the type of radar you are using. Is the radar for surface targets only (2-D, as in many marine radars), or for tracking airborne targets (3-D)? You will need to specify the mean clutter reflectivity (an area, for 2D, or a volume, for 3D). If it is a Doppler radar, you need to specify the clutter's Doppler spectrum. The radar signal processing algorithms greatly influence the clutter that is actually observed.

A very simple model, for a 2D non-Doppler radar:

Clutter returns are independent of one another within each sweep, and independent from one sweep to the next. The number of clutter returns per unit area, on each sweep, has a Poisson distribution, with mean λ. The locations of clutter returns are given by a 2D uniform distribution. The amplitude of clutter returns has an exponential distribution, with mean area .

Implementation of the very simple model:

lambda=10;                          % mean clutter returns per unit area [km^-2]
sigma0=1e-3;                        % mean clutter return area [km^2]
range=1;                            % radar range [km]
N1=poissrnd(lambda)*(2*range)^2;    % number of clutter returns in 2*range square
x1=-range+2*range*rand(1,N1);       % clutter x-coords [km]
y1=-range+2*range*rand(1,N1);       % clutter y-coords [km]
a1=exprnd(sigma0,1,N1);             % clutter areas [km^2]
r1=(a1/pi).^.5;                     % clutter radii [km]
x=x1(x1.^2+y1.^2<=range);   % clutter returns in range (exclude far corners)
y=y1(x1.^2+y1.^2<=range);   % clutter returns in range (exclude far corners)
r=r1(x1.^2+y1.^2<=range);   % clutter returns in range (exclude far corners)
N=length(x);                % number of clutter returns in range
%% Plot clutter returns
figure
for i=1:N
    rectangle('Position',[x(i)-r(i),y(i)-r(i),2*r(i),2*r(i)],...
              'Curvature',[1,1],'Facecolor','b','EdgeColor','none')
end
% Add range circles
rectangle('Position',0.25*[-1,-1,2,2],'Curvature',[1,1],'EdgeColor','r')
rectangle('Position',0.50*[-1,-1,2,2],'Curvature',[1,1],'EdgeColor','r')
rectangle('Position',0.75*[-1,-1,2,2],'Curvature',[1,1],'EdgeColor','r')
rectangle('Position',1.00*[-1,-1,2,2],'Curvature',[1,1],'EdgeColor','r')
xlim([-range,range]); ylim([-range,range])
axis equal

The clutter amplitude model above (which is that the clutter return areas have an exponential distribution) does not allow us to control the mean and variance independently. A gamma distribution (a scaled chi-squared distribution) allows you to set the mean and the variance independently.

Good luck!

[Kerem Yüksel]

At first, thanks a lot for your effort. How can I analyze the signal of that clutter model?

Have a good day, thanks.

[William Rose]

@Kerem Yüksel,

You can't. Which is why I said it is a very simple model.

Are you not using thre radar toolbox? If you are not using the toolbox, then I assume you have written your own code to simulate and analyze radar returns. The modeling of clutter will depend on the nature and level of your simulation of the radar return signal. I recommend the book I mentioned above, and Radar Signal Analysis and Procressing in Matlab, which does not appear to use the Radar Toolbox. Available online here. See chapter 9, . Section 9.5 and equations 9.45-9.48, and 9.55 and following. seem especially relevant.

[Kerem Yüksel]

assume you have written your own code to simulate and analyze radar returns. -> no I don' t but looking to write.

[William Rose]

@Kerem Yüksel,

I would first write code to simulate and process radar returns without clutter. I would test my code to become confident that it is working as expected. Then I would consider adding clutter to the simulation. I expect this will be a substantial effort. I hope the two books I have recommended above will be helpful.

[Subhosri]

@William Rose I am trying to model heterogenous clutter from airborne, by taking random samples from K,Weibull, lognormal and gamma distribution and I am trying to get a Range-Doppler map. But my clutter signal is always strongly concentrated at Doppler 0. Can you help me in figuring out what I am doing wrong?

d = 0.5;          % Array spacing [m]
va = 150;         % Platform speed [m/s]
phi_p = -pi/2;    % Array axis vs motion direction (forward-looking)
N = 4;            % Array elements
Na = 2;           % Range ambiguity factor
Nc = 2;         % Clutter scatterers
l = 1;            % Range bin index
f= 1*10^9;       %  L Band Radar Frequency
R = (10000*rand(5,1)+2000);% Range
u = 100 * ( rand(5,1) - 0.5 ); % Speed
c=3*10^8;%Speed of Light
Fd= ((2*u)/c)*f;% Doppler Frequency 
% Slow Time Bins
M = 1024; %no.of pulses
L = 1000;% Fast Time Bins
PRF=1000;% Slow Time PRF 
Tm = 1/PRF;% Slow Time Period
Fs=10*10^6;% Fast Time Sampling rate 
tl= 1/Fs;% Fast time period
p = [ 1 0.5 0.25 0.25 0.1 ]';% Target reflectivity
barker = [1 1 1 1 1 -1 -1 1 1 -1 1 -1 1];% Barker code for 13 chips
Ko = ( (2*pi)/.30 );% Propagation Constant
% Filter Coefficients % conjugate
a = [1 -1 1 -1 1 1 -1 -1 1 1 1 1 1 ];
% Unambiguous range for doppler 
Run = c/(2*PRF);
% Time Delay Ranging for Monostatic Radar
to= round ( 2*R/(3*10^8)/tl);
% Signal to clutter ratio
SCR = 100;
%clutter modeling
C = simulate_clutter(M, N, Na, Nc, l, f, c, PRF, d, va, phi_p).';

for i = 1 : M
    % The Time Dependence
    t=i*Tm;
    % Copying Clutter to the Signal
    S=C;
    % The Echo received for each Target +  Clutter
    for k = 1 : 5
        Delay = to(k);
        S(Delay:Delay+12) = S(Delay:Delay+12)+ p(k).* barker (1:13).* exp(j*2*Ko*u(k)*t);
    end
    %Matched Filtering with inverse Barker code coefficients
    Sfilterd= filter(a,1,S)/13;
    % Building the Matrix Cube
    MatrixCube(i,1:L) = Sfilterd(1:L);
end
FFTMatrixCube = fftshift(fft(MatrixCube,[],1));
figure 
% Ploting the Doppler - Range Grapgh
image([-500 500],[-15 15],200*log10(abs(FFTMatrixCube')));
xlabel('Doppler')
ylabel('Range')


function xC = simulate_clutter(K, N, Na, Nc, l, fc, c, fr, d, va, phi_p)
    T = 1/fr;
    xC = zeros(K*N, 1);  %clutter signal

    dist_list = {'k', 'gamma', 'weibull', 'lognormal'};  % Distribution pool
    for m = 1:Na
        for n = 1:Nc
            phi = rand() * 2 * pi;       % azimuth
            varphi = rand() * pi/2;      % elevation

            cos_psi = cos(phi - (pi + phi_p)) * cos(varphi);
            cos_theta = cos(phi) * cos(varphi);

            f_dm = 2 * va * cos_psi * fc / (c * fr);     % Doppler freq
            f_sm = cos_theta * fc * d / c;               % spatial freq

            s_t = exp(1j * 2 * pi * f_dm * (0:K-1)).';    % K x 1
            s_s = exp(1j * 2 * pi * f_sm * (0:N-1)).';    % N x 1

            dist_choice = dist_list{randi(length(dist_list))}; 
            xi = generate_random_xi(dist_choice, K * N);
            xi = apply_temporal_correlation(xi);  % K*N x 1

            xC = xC + kron(s_t, s_s).*xi;
        end
    end
end

function xi = generate_random_xi(type, K)
    CNR = 60 - 10.*log10(4*1024);
    Pc = 10^(CNR/10); 
    switch lower(type)
        case 'k'
            g = randraw('gamma', [1/3], K);
            r = randraw('rayl', 1, K);
            xi = sqrt(g) .* r;
            xi = xi./sqrt(xi'*xi).*sqrt(Pc); 
        case 'weibull'
            xi = randraw('weibull', [0, 1, 2], K);
            xi = xi./sqrt(xi'*xi).*sqrt(Pc); 
        case 'gamma'
            xi = randraw('gamma', [1/3], K);
            xi = xi./sqrt(xi'*xi).*sqrt(Pc); 
        case 'lognormal'
            xi = randraw('lognorm', [-1, 10.2], K );
            xi = xi./sqrt(xi'*xi).*sqrt(Pc); 
        otherwise
            error('Unknown distribution type.');
    end
end

function xi_corr = apply_temporal_correlation(xi)
    rho = 0.5; 
    K = length(xi);
    xi_corr = zeros(size(xi));
    xi_corr(1) = xi(1);
    for k = 2:K
        xi_corr(k) = rho * xi_corr(k-1) + sqrt(1 - rho^2) * xi(k);
    end
end

[William Rose]

@Subhosri,

Clutter returns that are strongly clustered around Doppler 0 is what I would expect for the simple model of clutter returns which I described in my previous post, because those clutter returns are from non-moving objects.

There are many aspects of your code which I do not understand. These aspects are probably correct, but I do not know enough about Doppler radar processing to fully assess your code, and I do not have time to learn. Therefore my comments below will not solve your problem, but I hope they are useful.

I you have not done so already, I recommend that you test your main program by verifying that it gives the expected results when there are no clutter returns. Do this by adding a line at the end of function simulate_clutter(), which will set xC to be a vector of zeros. Comment out this line after you have verified that your main program is otherwise working as you expect.

It appears that Nc=number of clutter scatterers=2. Is that correct? Does it mean there are 2 clutter sources? If so, that is not what I would expect. I would expect the number of clutter retuns to be a discrete random variable, perhaps with a Poisson distribution.

It appears that simulate_clutter() choose a random location for each of the Nc clutter returns by seeting phi (azimuth) and varphi (elevation). YOur elevation random variable is uniform on [0,pi/2]. This means the clutter location is lust as likely to be between 89 and 90 degrees elevation as it is to be between 0 and 1 degree elevation. This is unrelaistic, because the area of the 0-to-1 degree band around a hemisphere is much larger than the area of the 89-to-90 degree disk at the top of a hemisphere. If you visualize the reuslting points with N=100 to 500 points, you will see that they are overly clustered around the north pole. A better way is to find a random point in the [-1,+1] cube, then texst to see it it is inside the unit sphere. If it is not, discard it. If it is inside, then find its azimuth and elevation, and that is your random point.

simulate_clutter() sets a Doppler frequency using va=platform velocity. Although I don ot understand all the details of your code, I suspect that this means the Doppler frequency of the clutter return will be what is expected for a stationary object. This may be why your clutter returns are strongly clustered arounf doppler=0. If you wan the clutter returns to have random velocities, then define a velocity range or standard deviation, then choose a random number for the clutter velocity. You may need to add this to va. Then use this value to calculate the Doppler frequency of the clutter return.

Analysis of the code:

The model of clutter returns is in function simulate_clutter():

function xC = simulate_clutter(K, N, Na, Nc, l, fc, c, fr, d, va, phi_p) 

This function returns a clutter return waveform, a K*N-by-1 vector. The return from actual targets will be added to the clutter return waveform in the main program. xC is initilaized to zero in the function. xC is updated NaxNc times, inside two nested for loops, where Na=range ambiguity factor and Nc=number of clutter scatterers. On each pass, a new vector is added to the existing xC vector. The new vector that is added on each pass is computed by calling generate_random_xi() and apply_temporal_correlation().

I encourage you to document your functions by including comments that deifne the meaning and size of all inputs and outputs. This will make it easier for you to understand your own code (especially if you come back to it after several months), find mistakes, and it it very important if you want to get help from other people, so they do not have to figure out what each variable is. In this case, the documentation for simulate_clutter() would look like this:

function xC = simulate_clutter(K, N, Na, Nc, l, fc, c, fr, d, va, phi_p) 
% Generate clutter return waveform. The return from targets will be added
% to this waveform.
% Inputs
%  K=number of pulses
%  N=number of array elements
%  Na=range ambiguity factor 
%  Nc=number of clutter scatterers
%  l=range bin index (scalar)
%  fc=L band radar frequency [Hz]
%  c=speed of light [m/s]
%  etc...
% Output
%  xC=K*Nx1 array of double [arbitrary units]

I also strrongly encourage you to include units in your main progeam every time you define a constant or variable that has units - including times, frequencies, velocities, and distances.

Good luck with your work.