%——無調頻 矩形 相參脈沖串 的模糊函數 called——

%模糊函數公式參考源自：書《雷達原理(第一版)》P1-101的式(140)和P1-95的式(109)。
function x = rectangle_pulses_ambig_called(taup,Tr,N,fdy,taux)

nf = length(fdy);
nt = length(taux);
x = zeros(nf,nt);%給模糊函數所有點賦初值0

fdmin = fdy(1);
fdmax = fdy(nf);
delta_fd = fdy(2) - fdy(1);

taumin = taux(1);
taumax = taux(nt);
delta_tau = taux(2) - taux(1);

jj = 0;
for tau = taumin : delta_tau : taumax
    jj = jj + 1;
    ii = 0;
    for fd = fdmin : delta_fd : fdmax
        ii = ii + 1;
        for p = -(N-1) : 1 : (N-1)
            val1 = exp(j * p * pi * fd * Tr);
            if (abs(tau - p * Tr) <= taup)
                t_temp = tau - p * Tr;
                x0 = exp(j * pi * fd * t_temp) * ((taup - abs(t_temp)) / taup) * sinc(fd * (taup - abs(t_temp)));
            else
                x0 = 0;
            end
            if  fd ~= 0
                val2 = sin(pi * fd * Tr * (N - abs(p))) / (N * sin(pi * fd * Tr));
            else
                val2 = (N -abs(p)) / N;
            end
            temp(p + (N - 1) + 1) = val1 * x0 * val2;
        end
        x(ii,jj) = abs(sum(temp));
    end
end


%公式三
% function x = rectangle_pulses_ambig_called(taup,Tr,N,taumin,delta_tau,taumax,fdmin,delta_fd,fdmax)
% eps = 0.000000001;
% dimension_i = (fdmax - fdmin) / delta_fd + 1;
% dimension_j = (taumax - taumin) / delta_tau + 1;
% x = zeros(dimension_i,dimension_j);%給模糊函數所有點賦初值0
% for p = -(N-1) : 1 : N-1
%     j = 0;
%     for tau = taumin :delta_tau : taumax
%         j = j + 1;
%         if (tau >= p*Tr-taup) & (tau <= p*Tr+taup)
%             i = 0;
%             for fd = fdmin : delta_fd : fdmax
%                 i = i + 1;
%                 val1 = pi * fd * Tr * (N - abs(p));
%                 val2 = pi * fd * (taup - abs(tau - p * Tr));
%                 val3 = pi * fd * Tr;
%                 val4 = N * pi * fd * taup;
%                 val5 = exp(j * pi * fd * tau);
%                 x(i,j) = x(i,j) + val5 * (sin(val1) * sin(val2) + eps) / (sin(val3) * val4 + eps);
%             end
%         end
%     end
% end
% x=abs(x);

%公式二
% function x = rectangle_pulses_ambig_called(taup,Tr,N,taumin,delta_tau,taumax,fdmin,delta_fd,fdmax)
% eps = 0.000000001;
% j = 0;f(1:2,1:3)=0
% for tau = taumin : delta_tau : taumax
%     j = j + 1;%外循環變量為tau，模糊函數x(i,j)應沿著x軸變化,所以外循環中的增量應為j
%     i = 0;
%     if abs(tau) < taup
%         for fd = fdmin : delta_fd : fdmax
%             i = i + 1;%內循環變量為fd，模糊函數x(i,j)應沿著y軸變化,所以內循環中的增量應為i
%             acc = 0;
%             for p = -(N-1) : 1 : N-1
%                 val1 = sin(pi * fd * Tr * (N - abs(p)) + eps) / (N * sin(pi * fd * Tr) + eps);
%                 val2 = taup - abs(tau + p * Tr);
%                 val3 = val2 / taup;
%                 val4 = sin(pi * fd * val2 + eps) / (pi * fd * val2 + eps);
%                 acc = acc + val1 * val3 * val4;        
%             end
%             x(i,j) = abs(acc);
%             endfgrg  
%     end
% end

%公式一
% function x = rectangle_pulses_ambig_called(taup,Tr,N,taumin,delta_tau,taumax,fdmin,delta_fd,fdmax)
% eps = 0.00000000001;
% j = 0;
% for tau = taumin : delta_tau : taumax
%     j = j + 1;%外循環變量為tau，模糊函數x(i,j)應沿著x軸變化,所以外循環中的增量應為j
%     i = 0;
%     for fd = fdmin : delta_fd : fdmax
%         i = i + 1;%內循環變量為fd，模糊函數x(i,j)應沿著y軸變化,所以內循環中的增量應為i
%         acc = 0;
%         for p = -(N-1) : 1 : N-1
%             val1 = pi * fd * (taup - abs(tau - p * Tr));
%             val2 = pi * fd * (N - abs(p)) * Tr;
%             val3 = N * taup * pi * (fd + eps) * sin(pi * (fd + eps) * Tr);
%             val4 = exp(j * pi * fd * (N - 1 + p) * Tr) *  exp(j * 2 * pi * fd * (tau - p * Tr));
%             acc = val4 * sin(val1) * sin(val2) / val3 + acc;
%         end
%         x(i,j) = abs(acc);
%     end
% end