helpful.
I think the code has little bug. It does not work if we set doppler shift (fd) is very small.
I made a small correction, and it works fine for me. just add one line.
fM=fd/fs; %normalized doppler shift
NfM=fftSize*fM;
kM=floor(NfM); %maximum freq of doppler filter in FFT samples
if (kM < 1) kM = 1; end % newly added line
% rest of the code
% ...
I am not sure though why it is like this, as I am not expert in channel analysis.