% 雜波散射點(diǎn)的仿真
clear all;clc;

% 基本參數(shù)設(shè)定
K = 3;         % 脈沖的個(gè)數(shù)
M = 4;         % 雜波散射點(diǎn)個(gè)數(shù)
% Texture 產(chǎn)生
b = 1;   % scale 參數(shù)
v = 0.5; % shape 參數(shù)
T = gamrnd(v,b,M,1);   % 對(duì)應(yīng) M 個(gè)散射點(diǎn)的Texture
% speckle 產(chǎn)生
a = abs(normrnd(0,0.1,K,M));    % 雜波散射點(diǎn)的反射率 
b = abs(randn(1,K)');      % 目標(biāo)對(duì)入射波的反射率
tau_t = 3 ;           % 目標(biāo)的延遲時(shí)間
tau_c = [1,5,7,8];    % 雜波散射點(diǎn)的延遲時(shí)間 
% 發(fā)射信號(hào)
fs = 100;       % 采樣頻率  hz
ts = 1/fs;       % 采樣間隔  us
np1 = 0.3;
np2 = 0.6;
taup = np2 - np1 ;  % 脈沖寬度
ns = 0;
nf = 1;
[x1,n] = stepseq(np1,ns,nf,ts);
[x2,n] = stepseq(np2,ns,nf,ts);
s = x1 - x2;        % 發(fā)射信號(hào)
ssp = zeros(K,length(n));
nssp = zeros(K,length(n));
prf = 10;            % 脈沖重頻
ss = s;
nss = n;
for i = 1:K
    [ssp(i,:),nssp(i,:)] = seqshift(s,n,prf*i);
    ss = [ss,ssp(i,:)];
    nss = [nss,nssp(i,:)];
end
plot(nss,ss);        % 最終發(fā)射信號(hào)

% 回波信號(hào)
[y_t,ny_t] = seqshift(ss,nss,tau_t);   % 目標(biāo)回波信號(hào)
nlegth = length(n);
for j =1:K
    y_t(((j-1)*nlegth+1):(j*nlegth)) = y_t(((j-1)*nlegth+1):(j*nlegth))*b(j,1);
end
y_c = zeros(M,length(ss));             % 雜波散射點(diǎn)回波信號(hào)
ny_c = zeros(M,length(ss));
for i = 1:M
    [y_c(i,:),ny_c(i,:)] = seqshift(ss,nss,tau_c(i));
    for j =1:K
        y_c(i,((j-1)* nlegth+1):(j*nlegth)) = y_c(i,((j-1)*nlegth+1):(j*nlegth)) * a(j,i);
    end
end
nr = zeros(K,length(nss));
sr = zeros(K,length(nss));        % 回波信號(hào)
% 構(gòu)造回波數(shù)據(jù)矩陣
for j = 1:K
    ons = y_c(:,((j-1)*nlegth+1):(j*nlegth));   % 回波數(shù)據(jù)
    nns = ny_c(:,((j-1)*nlegth+1):(j*nlegth));  % 回波位置信息
    onns = ons(1,:);
    nnns = nns(1,:);
    for i = 2:M
        onns = [onns,ons(i,:)];
        nnns = [nnns,nns(i,:)];
    end
    for p = 1:length(nnns)
        for q = (p+1):length(nnns)
            if nnns(p) == nnns(q)
                onns(p) = onns(p) + onns(q);
                nnns(q) = -100;
            end
        end
    end
    nrr = nnns(find(nnns~=-100));
    gr = onns(find(nnns~=-100));
end

% 
% nrp = [ny_c(1,:),ny_t,ny_c(2,:),ny_c(3,:),ny_c(4,:)];
% grp = [y_c(1,:),y_t,y_c(2,:),y_c(3,:),y_c(4,:)];
% for i=1:length(nrp)
%     for j =(i+1):length(nrp)
%         if nrp(i) == nrp(j)
%             grp(i) = grp(i) + grp(j);
%             nrp(j) = -100;
%         end
%     end
% end
% nr = nrp(find(nrp~=-100));
% gr = grp(find(nrp~=-100));
figure;plot(nr,gr);
zs = xcorr(ss,gr);
R = corrcoef(zs');
rrank = rank(R)