配列によって信号を保存するBPSK変復調シミュレーションにグラフを追加

list4_1.c

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

#define PI acos(-1) //3.1459
#define L 100//number of transmit bits
#define SNR_STEP 10

typedef struct
{
    double I;
    double Q;
}COMPLEX;
/*一様乱数生成*/
double _randU(void){
    return((double)rand()/(double)RAND_MAX);
}

double _randN(void){
    double s, r, t;
    s = _randU();
    if(s == 0.0) s = 0.000000001;
    r = sqrt(-2.0*log(s));
    t = 2.0*PI*_randU();
    return(r*sin(t));
}
/*ランダムデータ発生*/
void _randData(unsigned int *d, unsigned int length){
    unsigned int i;
    double x;
    for(i=0; i<length;i++){
        x = _randU();
        if(x >= 0.5)
            *(d+i) = 1.0;
        else
            *(d+i) = 0.0;
    }
}
/*AWGN発生*/
void _awgn(COMPLEX *n, unsigned int length, double Pn){
    unsigned int i;
    for(i=0; i<length;i++){
        (n+1)->I = _randN() * sqrt(Pn/2);
        (n+1)->Q = _randN() * sqrt(Pn/2);
    }
}
/*雑音電力計算*/
double _SNRdB2noisePower(double c_dB){
    return pow(10,(-1)*c_dB/10);
}

/*BPSK変調器*/
double _bpskMod(COMPLEX *s, unsigned int *data, unsigned int length){
    unsigned int i;
    for(i=0;length<i;i++)
        (s+i)->I = (*(data+i)-0.5)*2.0;
}
/*ベクトル加算*/
void _vectorSum(COMPLEX *a, COMPLEX *b, COMPLEX *c, unsigned int length){
    unsigned int i;
    for(i=0;i<length;i++){
        (a+i)->I = (b+i)->I + (c+i)->I;
        (a+i)->Q = (b+i)->Q + (c+i)->Q;
    }
}
/*BPSK復調器*/
void _bpskDem(unsigned int *rData, COMPLEX *rs, unsigned int length){
    unsigned int i;
    for(i=0;i<length;i++){
        if((rs+i)->I>0)
            *(rData+i) = 1;
        else
            *(rData+i) = 0;
    }
}
/*BER計算*/
double _BER(unsigned int *data,unsigned int *rData, unsigned int length){
    unsigned int i, sum = 0;
    double BER;
    for(i=0;i<length;i++)
        sum += abs(*(data+i) - *(rData+i));
    BER = (double)sum/length;
    return(BER);
}
/*得られたBERのファイルへ書き出し*/
void _BERprint(unsigned int snr_step, double *SNRdB, double *BER){
    unsigned int i;
    FILE *fp;
    if((fp = fopen("BER.txt","w"))==NULL){
        printf("Error:Can not open your file \n");
        exit(1);
    }
    else{
        for(i=0;i<snr_step;i++){
            printf("%f [dB] BER ~ %f\n", *(SNRdB+i), *(BER+i));
            fprintf(fp,"%f %f\n",*(SNRdB+i),*(BER+i));
        }
        fclose(fp);
    }
}

int main(void){
    FILE *gp;
    unsigned int txData[L], rxData[L], i, x[L];
    COMPLEX txSymbol[L], noise[L], rxSymbol[L];
    double BER[SNR_STEP], SNRdB[SNR_STEP];
    for(i=0; i<L; i++) x[i] = i;
    for(i=0; i<SNR_STEP; i++) SNRdB[i] = (double)i;
    for(i=0; i<SNR_STEP; i++){
        srand(time(NULL));

        _randData(txData,L);

        _bpskMod(txSymbol,txData,L);

        _awgn(noise,L,_SNRdB2noisePower(SNRdB[i]));

        _vectorSum(rxSymbol,txSymbol,noise,L);

        _bpskDem(rxData,rxSymbol,L);

        BER[i]=_BER(txData,rxData,L);
    }
    //_BERprint(SNR_STEP, SNRdB, BER);

    gp = popen("gnuplot -persist","w");//パイプを開き，gnuplotを立ち上げる
    fprintf(gp, "set multiplot\n");//マルチプロットモード
    fprintf(gp,"set xrange [0:%d]\n",L);//範囲の指定
    fprintf(gp,"set yrange [0:1]\n");//範囲の指定
    fprintf(gp, "set xlabel \"data\"\n"); // ラベル表示
    fprintf(gp, "set ylabel \"y\"\n");
    /* 送信データ */
    fprintf(gp,"plot '-' with lines linetype 1\n");
    for(i=0;i<L;++i){
        fprintf(gp,"%d\t%d\n",x[i],txData[i]);//送信データをプロット
    }
    fprintf(gp,"e\n");
    /* 受信データ */
    fprintf(gp,"plot '-' with lines linetype 3\n");
    for(i=0;i<L;++i){
        fprintf(gp,"%d\t%d\n",x[i],rxData[i]);//送信データをプロット
    }
    fprintf(gp,"e\n");

    fprintf(gp,"set nomultiplot\n");// マルチプロットモード終了
    fprintf(gp, "exit\n");//gnuplotの終了
    fflush(gp);
    pclose(gp);//パイプを閉じる
    return 0;
}