/* Pembangkitan tabel metrik untuk keperluan analisis kode konvolusi * dengan asumsi kanal AWGN berbasis sistem modulasi PSK * * Diawali dengan mengevaluasi fungsi probabilitas normal dan selanjutnya * menghitung fungsi log-likelihood untuk setiap nilai simbol yang * kemungkinan akan muncul * * Copyright 1995 Phil Karn, KA9Q * Updated March 2014 (!) */ /* Simbol-simbol itu adalah suatu nilai biner-ofset antara 0 - 255 dengan * asumsi nilai 128 adalah nilai dimana kemungkinan "0" dan "1" pada * posisi "bingung", atau lazim disebut suatu kondisi "erased" atau simbol * nirmakna atau nirinformasi. */ #define GNUPLOT "gnuplot -persist" #define notdef #define OFFSET 128 #include #include #include #include /* Integral fungsi Normal dari -Inf ke x. Rentang: 0-1 */ #define normal(x) (0.5 + 0.5*erf((x)/M_SQRT2)) // M_SQRT2 = akar 2 /* Logarithm base 2 */ #define log2(x) (log(x)*M_LOG2E) // log2(x) = ln(x) * ln(2) /* Pembangkit metrik log-likelihood untuk 8-bit soft-quantized-channel * dengan asumsi AWGN dengan modulasi BPSK */ int main() { int mettab[2][256]; /* Metric table, [sent sym][rx symbol] */ int wmettab[2][256]; /* mettab yang dipakai oleh wsprd */ double pdftab[2][256]; // int amp=100; /* Signal amplitude, units */ int amp=50; /* Signal amplitude, units */ double noise; /* Relative noise voltage */ // double bias=0.5; /* Metric bias; 0 for viterbi, rate for sequential */ double bias=0.45; /* Metric bias; 0 for viterbi, rate for sequential */ int scale=10; /* Scale factor */ int s,bit; double metrics[2][256]; double p0,p1,metrics0,metrics1; double mebn0 = 3.0; // double mebn0 = 3; double mesn0; char gambar_pdf[50]; char gambar_prob[50]; char gambar_wprob[50]; FILE *fpdf,*fprob,*fwprob,*gp,*gs; #include "./metric_tables.c" strcpy(gambar_pdf,"./pdf.tmp"); strcpy(gambar_prob,"./prob.tmp"); strcpy(gambar_wprob,"./wprob.tmp"); if((fpdf = fopen(gambar_pdf,"w")) == NULL) { fprintf(stderr, "Cannot open data file '%s'....\n", gambar_pdf); exit(1); } if((fprob = fopen(gambar_prob,"w")) == NULL) { fprintf(stderr, "Cannot open data file '%s'....\n", gambar_prob); exit(1); } if((fwprob = fopen(gambar_wprob,"w")) == NULL) { fprintf(stderr, "Cannot open data file '%s'....\n", gambar_wprob); exit(1); } if((gp = popen(GNUPLOT,"w"))==NULL){ /* 'gp' is the pipe descriptor */ printf("Error opening pipe to GNU plot. Check if you have it! \n"); exit(0); } if((gs = popen(GNUPLOT,"w"))==NULL){ /* 'gs' is the pipe descriptor */ printf("Error opening pipe to GNU plot. Check if you have it! \n"); exit(0); } mesn0 = mebn0 + 10*log10(bias); /* metric table Es/N0 in dB */ noise = sqrt(0.5/pow(10.,mesn0/10.)); /* Zero mempunyai makna tersendiri, yakni semua probabilitas antara -takberhingga hingga 0 "dititipkan" pada harga ini. Artinya adalah bahwa semua harga pada ekor sebelah kiri akan dititipkan pada zero. */ /* Prob s=0 volt jika dikirim digit "1" */ /* P(s|1) */ p1 = normal(((0-OFFSET+0.5)/amp - 1)/noise); /* Prob s=0 volt jika dikirim digit "0" */ /* P(s|0) */ p0 = normal(((0-OFFSET+0.5)/amp + 1)/noise); metrics[0][0] = log2(2*p0/(p1+p0)) - bias; metrics[1][0] = log2(2*p1/(p1+p0)) - bias; pdftab[0][0]=p0; pdftab[1][0]=p1; // -bias untuk kasus p0 = p1 maka metrics[0][s]=metrics[1][s]= -bias ??? // double plus,minus; for(s=1;s<255;s++){ /* P(s|1), prob diterima s volt jika dikirim digit "1" */ p1 = normal(((s-OFFSET+0.5)/amp - 1)/noise) - normal(((s-OFFSET-0.5)/amp - 1)/noise); // konstanta -1 menunjukkan area PDF-1 /* P(s|0), prob diterima s volt jika dikirim digit 0 */ p0 = normal(((s-OFFSET+0.5)/amp + 1)/noise) - normal(((s-OFFSET-0.5)/amp + 1)/noise); // konstanta +1 menunjukkan area PDF-0 if(p0 == p1){ // metrics0 = metrics1 = -bias; /* wsprd apakah masih pakai ini ??? */ metrics0 = (p0 == 0) ? -33.0 : 1 + log2(p0) - log2(p1+p0)-bias ; metrics1 = (p1 == 0) ? -33.0 : 1 + log2(p1) - log2(p1+p0)-bias ; } else { metrics0 = (p0 == 0) ? -33.0 : log2(2*p0/(p1+p0)) - bias; metrics1 = (p1 == 0) ? -33.0 : log2(2*p1/(p1+p0)) - bias; // Ekivalen: // jika p0==0 maka metrics=-33; // jika p0=!0 maka metrics=1+log2(p0)-log2(p1+p0)-bias metrics0 = (p0 == 0) ? -33.0 : 1 + log2(p0) - log2(p1+p0) - bias; metrics1 = (p1 == 0) ? -33.0 : 1 + log2(p1) - log2(p1+p0) - bias; } /* inilah yang manjadi kunci sukses soft-decision convolution decoder */ metrics[0][s] = metrics0; /*prob sebagai digit 0 jika diterima s volt*/ metrics[1][s] = metrics1; /*prob sebagai digit 1 jika diterima s volt*/ pdftab[0][s]=p0; pdftab[1][s]=p1; } /* 255 juga mempunyai nilai khusus */ /* P(s|1) */ p1 = 1 - normal(((255-OFFSET-0.5)/amp - 1)/noise); /* P(s|0) */ p0 = 1 - normal(((255-OFFSET-0.5)/amp + 1)/noise); // printf("p(255|1)= %f; p(255|0)= %f\n",kanan1,kanan0); // printf("255 %f %f\n",kanan1,kanan0); metrics[0][255] = log2(2*p0/(p1+p0)) - bias; metrics[1][255] = log2(2*p1/(p1+p0)) - bias; pdftab[0][255]=p0; pdftab[1][255]=p1; #ifdef notdef /*Probabilitas eror (Pe) dari suatu simbol (level tegangan s) adalah probabilitas sebagai digit "1" untuk level tegangan s antara 0-128. Ini adalah sama dengan hasil integrasi dari kurva normal-ofset dari -inf ke 0. */ //printf("simbol Pe = %lg\n",normal(-1/noise)); // getchar(); #endif for(bit=0;bit<2;bit++){ for(s=0;s<256;s++){ /* diskala mendekati nilai integer terdekat */ mettab[bit][s] = floor(metrics[bit][s] * scale + bias); } } for(s=0;s<256;s++){ fprintf(fpdf,"%d %lf %lf\n",s,pdftab[0][s],pdftab[1][s]); } for(s=0;s<256;s++){ fprintf(fprob,"%d %d %d\n",s,mettab[0][s],mettab[1][s]); } fclose(fpdf); fclose(fprob); for(s=0; s<256; s++) { wmettab[0][s]=round( 10*(metric_tables[2][s]-bias) ); wmettab[1][s]=round( 10*(metric_tables[2][255-s]-bias) ); fprintf(fwprob,"%d %d %d\n", s,wmettab[0][s],wmettab[1][s]); } fclose(fwprob); fprintf(gp,"plot './pdf.tmp' using 1:2 lw 3 lt 1 with lines title 'P(s|0)'\n"); fflush(gp); fprintf(gp,"replot './pdf.tmp' using 1:3 lw 3 lt 3 with lines title 'P(s|1)'\n"); fprintf(gs,"plot './prob.tmp' using 1:2 lw 3 lt 1 with lines title 'P(0|s)'\n"); fflush(gs); fprintf(gs,"replot './prob.tmp' using 1:3 lw 3 lt 3 with lines title 'P(1|s)\n"); fflush(gs); fprintf(gs,"replot './wprob.tmp' using 1:2 lw 3 lt 7 with lines title 'WSPR P(0|s)\n"); fflush(gs); fprintf(gs,"replot './wprob.tmp' using 1:3 lw 3 lt 8 with lines title 'WSPR P(1|s)\n"); fclose(gp); fclose(gs); return 0; }