www.pudn.com > convolutive.rar > cps.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
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;