[C言語]ガウス性雑音をボックス・ミュラー法を用いて簡単に発生させるプログラミング[コード付き]

今回は確率信号を扱います。
まあガウス性雑音です。

全体のコードは一番下に貼ってあるので参考にどうぞ。

ガウス性雑音とは
ボックス・ミュラー法とは
コードの実装
グラフに出力してみた
全体のコード

ガウス性雑音とは

まあ簡潔に言うとノイズですね。

ノイズというと処理対象以外の不必要な情報のことやんな

その中でも、上にも書いてあるとおり取る値がガウス分布に乗っているということです。
つまり完全に適当に発生させたノイズではないってこと

今回はこのノイズを確率信号としてC言語で発生してみます。
その際にボックス・ミュラー法を用います。

ボックス・ミュラー法とは

まずガウス性雑音は、値xの発生する確率密度関数が下のf(x)で与えられる乱数です。

まあただのガウス分布の式やな

しかしボックス・ミュラー法を用いれば
以下の式のZで、ガウス性雑音が得られます。

おお　割とシンプルかも。
Zが二つあるけどこれはどういう意味？
あと分散と平均が式に入ってないけど大丈夫？

Zが二つあるのは1組の(X,Y)から２種類のガウス性雑音を発生できるということです。
つまり二つでるけど、どちらかを使っても、どっちもを使ってもおkってこと。

また分散と平均は式に入ってない分、式の段階では考慮できてないので、のちのコードの実装ではそれも考えてコードを書きます。

コードの実装

せっかく雑音を発生させたならグラフで見てみたいですよね。

なので今回は配列に入れた雑音をファイルに格納する関数も作ります。
ファイルに入れれば後にgnuplotなどでグラフを確認できますね。

ファイルに完成した雑音を格納する関数

まずこれから作りましょう。
Nは発生させる雑音の数で、マクロで定義しています。

とりあえずコード

void file(double G1[], double G2[]){

    FILE *fp;
    char f_name[64];
    printf("格納するファイル名を入力してください:");
    scanf("%s",f_name);
        
     if((fp=fopen(f_name,"w"))==NULL)
         printf("ファイルを開けませんでした\n");
     else{
         for(int i=0; i<N/2; i++){
             fprintf(fp,"%f\n",G1[i]);
             fprintf(fp,"%f\n",G2[i]);
         }
     }
    fclose(fp);
 }

仮引数で雑音が入った配列を受け取ってそれをファイルに格納すればよいですね。
仮引数が２つの配列になっているのは、上の式で書いたZ1とZ２両方受け取るためです。

よって格納する際のforループはN/2にしましょう。これで、配列が二つあるので格納される点はその倍のNになります。

ちなみにファイルについての記事を貼っておくので良かったらどうぞ。

雑音をrand関数をもちいて発生させる

まず上のボックス・ミュラー法の式をみて

logとsin,cosが出てくるのでmath.hをインクルード。
1/2乗はsqrt関数を用いる。
XとYはrand関数で発生させる。

っていう方針が立ちました。
やっていきましょう。

コードです

void BOX(double ave, double dis){
    double X,Y,Z1,Z2;
     double G1[N];
     double G2[N];

     for(int i=0; i<N/2; ){
         X = (double)rand()/RAND_MAX;
         printf("%f ",X);
         Y = (double)rand()/RAND_MAX;
         printf("%f\n",Y);
         
         if((X!=0)&&(Y!=0)){//log(0)でない場合のみノイズを計算する
            Z1 = sqrt(-2*log(X))*cos(2*M_PI*Y);
            Z2 = sqrt(-2*log(X))*cos(2*M_PI*Y);

            G1[i] =ave + dis*Z1;
            G2[i] =ave + dis*Z2;
            i++;//if文の中でだけインクリメント
         }
     }
    
    file(G1, G2);
}

まず一様乱数X,Yについては<stdlib.h>に、rand()関数という値をランダムで返す関数があるのでそれを利用します。

でもrand()関数で返ってくるランダム値の範囲って0~RAND_MAX(stdlib.hで定義されている)だから今回の0~1の範囲に全然あってないけど…

そうです。ただ普通に使うと0~1の範囲に収まらないので
rand()/RAND_MAXとすることで、値が0~1の間に収まります。
ちなみにRAND_MAXは処理系に依存しています。

またrand()もRAND_MAXもint型なので(double)でキャストしてます。

後はlogと三角関数、sqrt関数を使ってボックス・ミュラーの式をそのまま記述しました。
forループの回数はここでもN/2回です。

(追記2022/1/14)
rand関数ではXとYが０になる可能性がある、つまりlog(0)の計算を行う可能性がある指摘をいただきました。
なのでfor(int i=0; i<N/2; )として、
if文でXY共に0出ない時のみ処理とインクリメントを行い、
どう出ない時は何もせずにもう一度for文を回すことにしました。ご指摘ありがとうございました！！
コードはすでに訂正済みです。

計算した値を配列に格納する際に平均と分散を考慮しましょう。
平均と分散の値は仮引数で受け取ります。(上のコードだとaveとdis)
平均はその値に足して、分散はその値に掛けましょう。

それを配列に格納して上で作ったfile関数に渡しておｋです

グラフに出力してみた

ファイルに書き出した雑音をgnuplotでグラフにしてみました。

おｋ

てことで今回はここまで。
ではまた〜

全体のコード

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

#define N 64

void file(double G1[], double G2[]){

    FILE *fp;
    char f_name[64];
    printf("格納するファイル名を入力してください:");
    scanf("%s",f_name);
        
     if((fp=fopen(f_name,"w"))==NULL)
         printf("ファイルを開けませんでした\n");
     else{
         for(int i=0; i<N/2; i++){
             fprintf(fp,"%f\n",G1[i]);
             fprintf(fp,"%f\n",G2[i]);
         }
     }
    fclose(fp);
 }

void BOX(double ave, double dis){
    double X,Y,Z1,Z2;
     double G1[N];
     double G2[N];

     for(int i=0; i<N/2; ){
         X = (double)rand()/RAND_MAX;
         printf("%f ",X);
         Y = (double)rand()/RAND_MAX;
         printf("%f\n",Y);
         
         if((X!=0)&&(Y!=0)){//log(0)でない場合のみノイズを計算する
            Z1 = sqrt(-2*log(X))*cos(2*M_PI*Y);
            Z2 = sqrt(-2*log(X))*cos(2*M_PI*Y);

            G1[i] =ave + dis*Z1;
            G2[i] =ave + dis*Z2;
            i++;//if文の中でだけインクリメント
         }
     }
    
    file(G1, G2);
}

int main(int argc, const char * argv[]) {
    double ave,dis;
    
    printf("平均と分散を入力してください\n平均:");
    scanf("%lf",&ave);
    printf("分散:");
    scanf("%lf",&dis);
    
    BOX(ave, dis);

　　 return 0;
}