#include<stdlib.h>
#include<stdio.h>
#include<math.h>
#define pi 3.1415926535897931
#define SNR 21
#define M 2
FILE *file;
int N=int(2e6),dtc_s,er_s,i,j,k,dtc_i,dtc_q;
double m=double(RAND_MAX)/2.0;
double d,dic,I,Q,min,a,b,w,x,y,E_N,src_s;
double ps[SNR],map[M];
void main() {
for(i=0;i<SNR;i++) {
er_s=0;
E_N=pow(10,double(i)/10.);
d=sqrt(2*E_N/21);
//定義信號位置
map[0]=-1*d;
map[1]=d;
// 利用Uniform R.V.來決送那一個SYMBOL
for(j=0;j<N;j++) {
src_s=double(rand())/(RAND_MAX+0.01);
src_s=src_s*4;
// 產生GAUSSIAN R.V.
star: a=double(rand())/m-1;
b=double(rand())/m-1;
w=a*a+b*b;
if(w < 1 ) {
x=a*sqrt(-2*log(w)/w);
y=b*sqrt(-2*log(w)/w);
}
else goto star;
//不同symbol在I Q Channel產生不同的信號
// AWGN CHANNEL SIMULATION
I = map[int(src_s)%2] + x;
Q = map[int(src_s)/2] + y;
// detection Signal
//計算接受I Channel信號和各點之間的距離，用來判斷是那一個symbol。
min=100;
for(k=0;k<M;k++) {
dic = fabs(I - map[k]);
if(min > dic) {
min=dic;
dtc_i=k;
}
}
//計算接受到Q Channel信號和各點之間的距離，用來判斷是那一個symbol。
min=100;
for(k=0;k<M;k++) {
dic = fabs(Q - map[k]);
if(min > dic) {
min=dic;
dtc_q=k;
}
}
// I Q Chennel 的信號合成一個4-ary QAM symbol
dtc_s = 2*dtc_q+dtc_i;
if(dtc_s != int(src_s) )
er_s++;
}
// 計算 BER
printf("\n%d ",er_s);
ps[i]=double(er_s)/double(N);
}
//結果放入MATLAB的M-file中
file = fopen("C:\\MATLAB6p5\\work\\qam4.m","w");
fprintf(file,"a = [");
for(i=0;i<SNR;i++)
fprintf(file,"%.16f\n",ps[i]);
fprintf(file,"]';\n");
fprintf(file,"snr1=0:%d;\n",SNR-1);
fprintf(file,"semilogy(snr1,a,'k');hold on\n");
fclose(file);
}
