%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Rx = cps(x,T,hT,K,N,dx,verbose);
% calculate the cross-power-spectrum Rx

% check whether these parameters can be used with signals of length rx,
% this check is done by checking whether the last index that will be
% accessed during the calculation of Rx exists
if K*T*N > size(x,1)
error('the signals are too short for the chosen parameters')
end

% do the calculation
if verbose, disp('estimate the cross-power spectrum of x'), end
Rx = zeros(dx,dx,hT,K);
XX = zeros(dx,dx,hT);
for k=1:K
for n=1:N
if verbose, fprintf('%d/%d %d/%d \r',k,K,n,N); end
indx = (k-1)*T*N + (n-1)*T + (1:T);
X=fft(x(indx,:));
% since we are only interested in filter in frequency domain
% that correspond to real valued filters in time domain,
% we can ignore the negative frequencies
X=X(1:hT,:);
for row = 1:dx, for col = 1:dx
XX(row,col,:) = reshape(conj(X(:,col)).*X(:,row),[1 1 hT]);
end, end
Rx(:,:,:,k) = Rx(:,:,:,k) + XX;
end;
end;
Rx = Rx/N/T;