QPSK変調のマルチパスフェージングシミュレーション

前提・実現したいこと

MatlabによってQPSK変調のマルチパスフェージング（レイリーフェージング）下でのシミュレーションを実現したいです。
教科書ではマルチパス数（素波M）が8波程度でレイリー分布に従うようなのですがなかなかうまくいきませんどのようにすればうまくいくのか教えて頂きたいです。

該当のソースコード

MATLAB
1clear,close all
2% パラメータ設定
3Ndata = 10^6; %[bit](データビット数の設定)
4Nsymbol = Ndata/2 ; % シンボル数
5symrate = 5*10^6 ; %シンボルレート
6Ts = 1/symrate ; % シンボル周期
7Es_N0_dB = [0:20]; % Es/N0 
8% ランダムデータの生成
9dataSeq = MYrandData(Ndata);
10N = 8; %マルチパスの数
11
12c = 3.0*10^8; %[m/s]
13fc = 2.1*10^9; %[Hz]
14lambda = c/fc; %波長
15v = 40000/3600; %[m/s]
16fD = (v/lambda); %　最大ドップラー周波数
17
18% QPSK変調器
19
20spcOutput = reshape(dataSeq,2,Ndata/2); %(S/P変換)
21qpskSymbolIndex = [1,2]*spcOutput; %(10進数に変換)
22qpskSymbol = ones(1,Ndata/2) * exp(1j*pi/4);
23qpskSymbol(find(qpskSymbolIndex==1)) = exp(1j*3*pi/4);
24qpskSymbol(find(qpskSymbolIndex==3)) = exp(1j*5*pi/4);
25qpskSymbol(find(qpskSymbolIndex==2)) = exp(1j*7*pi/4);
26h_re = 0;
27h_im = 0;
28
29for ii = 1:length(Es_N0_dB)
30% AWGN通信路
31
32Pn(ii) = 10^(-(ii-1)/10); % 
33
34noise = (randn([1,Ndata/2]) + 1j * randn([1,Ndata/2])) *sqrt(Pn(ii)/2);
35
36for i = 1:N  
37% レイリーフェージング
38 theta = (2*pi)*rand(1,Ndata/2); %ランダム位相
39 phi = (2*pi)*rand(1,Ndata/2); %ランダム位相 
40%  ray = (1/sqrt(N))*(exp(1j*2*pi*fd+theta));
41h_re2 = (1/sqrt(N))*cos(2*pi*fD.*cos(phi)*Ts+theta); 
42h_im2 = (1/sqrt(N))*sin(2*pi*fD.*cos(phi)*Ts+theta);
43
44 h_re = h_re + h_re2 ;
45 h_im = h_im + h_im2 ;
46
47end
48
49ray = h_re + 1j*h_im ;
50%  r = real(qpskSymbol.*ray) + 1j*imag(qpskSymbol.*ray) + noise;
51r = qpskSymbol.*ray + noise;
52 
53r = (r.*conj(ray))./(abs(ray).^2);
54  
55%複調器
56demodData = zeros(2,Ndata/2);
57demodData(2,find(imag(r)<=0)) = 1;
58demodData(1,find(real(r)<=0)) = 1;
59qpskDemod = demodData(:);
60% BERの計算
61
62BER(ii) = sum(abs(dataSeq - qpskDemod))/Ndata;
63end
64
65% BERの理論値
66SNR = 10.^(Es_N0_dB/10);
67theoryBer = (1/2)*(1-sqrt(SNR./(SNR+1)));
68
69% 結果の表示
70figure
71semilogy(Es_N0_dB,theoryBer,'b.-','LineWidth',2);
72hold on
73semilogy(Es_N0_dB,BER,'mx-','Linewidth',2);
74% axis([0 20 10^-5 1])
75grid on
76legend('theory', 'simulation');
77% ylim([10^(-6) 1])
78xlabel('E/N, dB')
79ylabel('Bit Error Rate')
80title('Bit error probability curve for QPSK modulation') 
81