function H = channelpred(M, Nr, fs, Nfft, tau, pdp, fd, P, Dt, CPfrac, sigma2e)
R = (randn(Nfft*P,M,Nr)+j*randn(Nfft*P,M,Nr))/sqrt(2);
% P = repmat(sqrt(pdp),[1,M,Nr]);
% G = R.*P;
H = zeros(Nfft,M,Nr);
nidx = [-Nfft/2:Nfft/2-1]';
W = exp(-sqrt(-1)*2*pi*fs/Nfft*nidx*tau);
tsym = Nfft*(1+CPfrac)/fs;
rt = besselj(0,2*pi*fd*tsym*Dt*[0:P]);
Rt = toeplitz(rt(1:end-1));
Rh = W*diag(pdp)*W';
Rhat = kron(rt(2:end),Rh);
Rhathat = kron(Rt,Rh) + sigma2e*eye(P*Nfft);
[U,D] = eig(Rhathat);
Resqrt = Rhat*U*diag(sqrt(1./real(diag(D))));
for i = 1 : Nr
    H(:,:,i) = Resqrt*R(:,:,i);
end;