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

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

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



ガウス性雑音とは

ガウス雑音は正規分布と等しい確率密度関数を持つ統計的雑音。言い換えると、ノイズがとる値がガウス分布であるということである。 wikipedia[ガウス雑音]より引用


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

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



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

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

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

不思議ですね。
どうしてそうなるかはわからないですが、使わせてもらいましょう。

ただこれでは、分散と平均を考慮できてないので、のちのコードの実装ではそれも考えてコードを書きます。



コードの実装

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

なので今回は配列に入れた雑音をファイルに格納する関数も作ります。
ファイルに入れれば後に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);
 }


仮引数で雑音が入った配列を受け取ってそれをファイルに格納すればよいですね。

注意としては、
ボックス・ミュラー法では一組の一様乱数X,Yから二つのガウス性雑音Zが得られます。
なので受け取る配列は二つ。

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


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




雑音を発生させる関数

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

  • 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; i++){
         X = (double)rand()/RAND_MAX;
         Y = (double)rand()/RAND_MAX;

         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;
     }
    
    file(G1, G2);
}

まず一様乱数X,Yについて
0~1の範囲の乱数を発生させる必要があります。

<stdlib.h>に、rand関数という0~RAND_MAX(stdlib.hで定義されている)の範囲の値をランダムで返す関数があるのでそれを利用します。

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

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


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


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


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



グラフに出力してみた

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

おk


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

ではまた

#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; i++){
         X = (double)rand()/RAND_MAX;
         printf("%f ",X);
         Y = (double)rand()/RAND_MAX;
         printf("%f\n",Y);

         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;
     }
    
    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;
}

syumokuzame

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です

コメントする