通信路容量を導出するプログラム／ランダムな値が収束しない

Question

前提・実現したいこと:8PSKでの通信路容量の導出

Visual Studio 2019 でAWGM環境下で各変調ごとの通信路容量を導出するプログラムをC言語で書いています。BPSK,QPSKでの通信路容量は完璧に導出できたのですが、8PSKでの通信路容量のグラフの値がどれだけ多くの試行回数を重ねても収束しません。
BPSK,QPSKでは完璧にうまくいっているので、プログラムの内容自体はあっていると思うのですが、、、（つまり、私はこの問題は乱数が正規分布に収束していないという問題だと認識しています。）
どなたか助言いただける方がいたら、よろしくお願いします。
なお、エラーメッセージは一つも出ていません。
以下が、作りたいグラフです。

以下、出力時のグラフです。

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>
#include<complex.h>
#include"MT.h"
#pragma warning(disable : 4996)
#define modulation_level 8//BPSK:2 QPSK:4 8PSK:8 16QAM:16 変調多値数
#define N 10000000 //受信信号の要素数
#define EbN0db_max 20//プロットするEs/N0[db]の最大値[db]
#define EbN0db_min -10 //プロットするEs/N0[db]の最小値[db]
#define EbN0db_plot 31 //プロットする点数
#define Ns 1//1シンボル当たりのサンプル数
#define PI 3.1415926
#define amplitude 1 //雑音と信号の比が重要なので、振幅は1で固定

int genpulse()
{
	return genrand_int32() & 1;
}

int generate_constellation_diagram()//コンスタレーション上のプロットを表す数値を出力
{
	int pulse_number=0;

		for (int i = 0; i < log2(modulation_level); i++)
		{
			int pulse = genpulse();
			pulse_number+=pow(2, i)* pulse;
		}

	return pulse_number;
}

_Dcomplex generate_noise(long double N0)//N0を引数にして雑音(noise_I+j*noise_Q)を返す関数
{
	long double sigma = pow((long double)N0 / 2, 0.5);
	long double u1 = genrand_real3(); //(0,1)乱数を作成
	long double u2 = genrand_real3();
	long double noise_I = sigma* pow(-2 * log(u1), 0.5)* cos(2 * PI * u2);
	long double noise_Q = sigma* pow(-2 * log(u1), 0.5)* sin(2 * PI * u2);
	_Dcomplex noise = _Cbuild(noise_I, noise_Q);

	return noise;
}

long double calc_entropy_function(long double N0,int constellation_number)
{
	_Dcomplex constellation_point[modulation_level];
	if (modulation_level == 2)
	{
		constellation_point[0] = _Cbuild(amplitude, 0);
		constellation_point[1] = _Cbuild(-amplitude, 0);
	}

	else if (modulation_level == 4)
	{
		constellation_point[0] = _Cbuild(pow(2, -1 * 0.5) * amplitude, pow(2, -1 * 0.5) * amplitude);
		constellation_point[1] = _Cbuild(pow(2, -1 * 0.5) * amplitude, -pow(2, -1 * 0.5) * amplitude);
		constellation_point[2] = _Cbuild(-pow(2, -1 * 0.5) * amplitude, pow(2, -1 * 0.5) * amplitude);
		constellation_point[3] = _Cbuild(-pow(2, -1 * 0.5) * amplitude, -pow(2, -1 * 0.5) * amplitude);
	}

	else if (modulation_level == 8)
	{
		constellation_point[0] = _Cbuild(amplitude ,0);
		constellation_point[1] = _Cbuild(amplitude/sqrt(2), amplitude / sqrt(2));
		constellation_point[2] = _Cbuild(0,amplitude);
		constellation_point[3] = _Cbuild(amplitude / sqrt(2), -1 * amplitude / sqrt(2));
		constellation_point[4] = _Cbuild(0,-1*amplitude);
		constellation_point[5] = _Cbuild( -1*amplitude / sqrt(2), -1*amplitude / sqrt(2));
		constellation_point[6] = _Cbuild(0, -1*amplitude);
		constellation_point[7] = _Cbuild(amplitude/sqrt(2),-1*amplitude/sqrt(2));
	}
	_Dcomplex noise=generate_noise(N0);
	_Dcomplex	constellation_noise = _Cbuild(creal(noise)+creal(constellation_point[constellation_number]), cimag(noise)+cimag(constellation_point[constellation_number]));//得られる点のコンスタレーション上の座標
	long double tmp2=0;
	for (int i = 0; i < modulation_level; i++)//分子の合計をtmp2に返す
	{
		if (i == constellation_number)
		{
			long double tmp = 0;
			tmp = cabs(noise);
			tmp2 += exp(-1 * tmp * tmp / N0);

		}

		else
		{
			long double tmp = 0;
			tmp = cabs(_Cbuild(creal(constellation_noise) - creal(constellation_point[i]), cimag(constellation_noise) - cimag(constellation_point[i])));
			tmp2 += exp(-1 * tmp * tmp / N0);
		}

	}

	long double tmp=0;

	tmp = cabs(noise);

	return (log2(tmp2)-log2( exp(-1 * tmp * tmp / N0)));
}

int main()
{
	init_genrand((unsigned)time(NULL));
	long double channel_capacity=0;//C
	long double entropy_function=0;//H(x|y)
	int constellation_number;

	FILE* fp;
	
	if ((fp = fopen("d.txt","w")) == NULL)
	{
		printf("cannot open!");
		exit(1);
	}

	for (int i = 0; i < EbN0db_plot; i++)
	{

		long double EsN0db =i*((long double)EbN0db_max - EbN0db_min) / ((long double)EbN0db_plot -1)+ EbN0db_min;//-10~20まで31点プロット

		long double N0 = pow(10, -1 * EsN0db / 10);

		constellation_number=generate_constellation_diagram();//コンスタレーション上のプロットを表す整数を出力

		entropy_function = 0;//（nI+j*nQ)の整数とN0を引数にしてH(x|y)を求める関数
			for (int j = 0; j < N; j++)
			{
				entropy_function+= calc_entropy_function(N0, constellation_number);
			}
		entropy_function /= N;

		channel_capacity = log2(modulation_level) - entropy_function;//チャネル容量

		fprintf(fp, "%f,%lg\n", EsN0db, channel_capacity);

		printf("finish Es/N0=%f\n", EsN0db);
	}

	fclose(fp);

	return 0;

}

試したこと

・乱数生成が一様分布でない
→乱数生成の関数をrand関数からメルセンヌツイスタへ変更
・計算結果がdouble型に収まっていない
→すべてのdouble型をlong　double型に変更（complex.hに問題があるのかも？）
・測定回数が十分でない
→最大1000000000点で測定(終了まで15時間かかりました、、、)

補足情報（FW/ツールのバージョンなど）

実行環境:win10
visual studio 2019

Comments

[SHOMI]

>すべてのdouble型をlong　double型に変更
VisualC++のdoubleとlong doubleは同じサイズですよ。
https://docs.microsoft.com/en-us/cpp/cpp/data-type-ranges
long double same as double

[kaneko_]

ありがとうございます。もしサイズを増やしたかったら別の開発環境でコーディングしなきゃいけないですかね、、、

Answer 1

constellation_point[3]
constellation_point[4]

が間違っていました????
初歩的な間違いで申し訳ございませんでした????

Answer 2

calc_entropy_functionの8PSKのときの分岐で

    else if (modulation_level == 8)
    {
        constellation_point[0] = _Cbuild(amplitude ,0);
        constellation_point[1] = _Cbuild(amplitude/sqrt(2), amplitude / sqrt(2));
        constellation_point[2] = _Cbuild(0,amplitude);
        constellation_point[3] = _Cbuild(amplitude / sqrt(2), -1 * amplitude / sqrt(2));
        constellation_point[4] = _Cbuild(0,-1*amplitude);
        constellation_point[5] = _Cbuild( -1*amplitude / sqrt(2), -1*amplitude / sqrt(2));
        constellation_point[6] = _Cbuild(0, -1*amplitude);
        constellation_point[7] = _Cbuild(amplitude/sqrt(2),-1*amplitude/sqrt(2));
    }

となってますが、constellation_point[4]がおかしいような気がします。

constellation_point[4] = _Cbuild(-1*amplitude, 0);

じゃないですか?

Comments

[kaneko_]

本当ですね、間違っていました????
ごめんなさい
ですが、指摘箇所を変更してもやはり変わらず収束しません、、