Thanks a lot.
%%
Nt = 4; Nr = 4; % transmit and receive antennas
M = 16; L = 16; % CPI pulses number; code length
Fr = 1000; F0 = 1e9; B = 1e6; % PRF; carrier frenquency; bandwidth
H = 9000; Vp = 150; % heigh and velocity of the plane
Range = 12728; ThetaT = 0; % range and azimuth of the target
PhiT = asin(H/Range); Vt = 45; % elevation and speed of target
Nc = 361; % clutter patches
PhiC = PhiT; % elevation of clutter, equals to target
ThetaC = linspace(0,360,361)/180*pi; % azimuth of each clutter patch
SNR = 0; INR = 30; % signal-to-noise-ratio, signal-to-intereference-ratio
%%
c = 3e8; lamda = c / F0; % speed of light, wavelength
dt = lamda * 2; dr = lamda / 2; % antenna space
ut = (0:(Nt-1)).'*dt; ur = (0:(Nr-1)).'*dr; % positions of transmit and receive antenna
% =========================== target===================
Ft = 2 * Vt / (lamda * Fr); % doppler of target
uFt = exp(1i*2*pi*(0:(M-1)).'*Ft); % temporal steer vector
at = exp(1i*(2*pi/lamda)*ut*cos(PhiT)*sin(ThetaT)); % transmit steer vector
bt = exp(1i*(2*pi/lamda)*ur*cos(PhiT)*sin(ThetaT)); % rceiver steer vector
vt = kron(kron(uFt,bt),at); % virtual steer vector
THETA = 10^(SNR/10)*(vt*vt');
% ======================== intereference ====================
temp_vc = zeros(M*Nr*Nt,M*Nr*Nt);
for ii = 1:Nc
Fc_i = 2*Vp*cos(PhiC)*sin(ThetaC(ii))/(lamda*Fr);
uFc_i = exp(1i*2*pi*(0:(M-1)).'*Fc_i); % temporal steer vector
ac_i = exp(1i*(2*pi/lamda)*ut*cos(PhiC)*sin(ThetaC(ii))); % transmit steer vector
bc_i = exp(1i*(2*pi/lamda)*ur*cos(PhiC)*sin(ThetaC(ii))); % receiver steer vector
vc_i = kron(kron(uFc_i, bc_i), ac_i); % virtual steer vector
temp_vc = temp_vc + vc_i*vc_i';
end
PHI = 10^(INR/10)*temp_vc;
% =======================================================================
r = Nt*L; % dimension of signal
p = M*Nr*L; % dimension of filter
tuple.w = euclideancomplexfactory(p,1); % complex eculian product
tuple.s = complexcirclefactory(r,1); % complex circle product
Manifold = productmanifold(tuple); % product manifold
problem.M = Manifold; %
problem.cost = @cost; % cost function
problem.egrad = @grad; % gradient
options.maxiter = 100;
options.beta_type = 'P-R';
options.minstepsize = 1e-6;
options.linesearch = @linesearch;
% checkgradient(problem);
[X, Xcost, info] = conjugategradient(problem,[],options);
% [X, Xcost, info] = steepestdescent(problem,[],options);
figure(1)
subplot(1,2,1)
plot([info.gradnorm],'-*')
ylabel('norm of grad')
xlabel('iter num')
title('curve of gradient')
subplot(1,2,2)
plot(10*log10(abs(-[info.cost])))
ylabel('cost function')
xlabel('iter num')
title('curve of cost function')
%% ==================== functions =====================
function f = cost(X) % equation(8)
w = X.w; s = X.s;
W = reshape(w,[L,M*Nr]);
Wwave = kron(W.', eye(Nt)); % W~
f1 = s'*Wwave'*THETA*Wwave*s; %
f2 = s'*Wwave'*PHI*Wwave*s + w'*w; %
f = - real(f1 / f2);
end
function g = grad(X) % ecudian gradient
w = X.w; s = X.s;
W = reshape(w,[L,M*Nr]);
Wwave = kron(W.', eye(Nt)); % W~
Sstar = reshape(s,[Nt, L]); % S*
Swave = kron(eye(Nr*M),Sstar); % S~
P = Swave'*THETA*Swave; Q = Swave'*PHI*Swave;
k = w'*Q*w + w'*w;
R = Wwave'*THETA*Wwave; U = Wwave'*PHI*Wwave;
belta = s'*U*s + w'*w;
egradw = -(2/k)*P*w + (2/k^2)*w'*P*w*(Q*w+w); % equation(22)
egrads = -(2/belta)*R*s + 2/(belta^2)*s'*R*s*(U*s); % equation(23)
egrad.w = egradw;
egrad.s = egrads;
g = Manifold.egrad2rgrad(X, egrad);
end
end