function H = channel(M, Nr, fs, Nfft, tau, pdp)
R = (randn(length(pdp),M,Nr)+j*randn(length(pdp),M,Nr))/sqrt(2);
P = repmat(sqrt(pdp),[1,M,Nr]);
H = zeros(Nfft,M,Nr);
nidx = [-Nfft/2:Nfft/2-1]';
W = exp(-sqrt(-1)*2*pi*fs/Nfft*nidx*tau);
for i = 1:Nr,
    H(:,:,i) = W*(R.*P);
end;