BPSK，QPSK的C语言仿真

本文首先给出BPSK的原理，再给出对应的C语言设计过程，并附代码。
当给出BPSK的过程后，QPSK即为分路问题，将直接给出代码。

一、什么是BPSK

BPSK可以被称为二进制相移键控，其利用载波的相位变化来传递数字信息，而振幅和频率保持不变。以初始相位0和π\piπ来表示二进制比特（符号）1和0。在传输时可以用以下波形函数表示：

s(t)=Asin⁡(ωt+φn)s\left( t \right) = A\sin \left( {\omega t + {\varphi _n}} \right)s(t)=Asin(ωt+φn)

其中，φn\varphi _nφn表示第nnn个符号的绝对相位，取值规则与之前所述相同。

二、BPSK仿真的问题处理

1. 流程确认

首先确定数据流的传递方式。其流程可以由以下过程表示：

由以上图片可知，仿真的难点在于信噪比和高斯噪声的加入。其他方式由低通等效原理可以省略高频调制步骤，至于其余步骤均为映射和判断的if else语句没有难度。

2. Ｃ语言随机信源的产生方式

首先应该了解C语言的随机数产生方式。
在标准C语言库stdlib.h内有rand()函数用于生成一个介于0与RAND_MAX之间的int型整数。其中RAND_MAX是一个整数，它与系统有关。
rand函数没有输入参数，由于rand函数在产生随机数前，需要系统提供的生成伪随机数序列的种子，rand根据这个种子的值产生一系列随机数。如果系统提供的种子没有变化，每次调用rand函数生成的伪随机数序列都是一样的，所以说C语言的随机是一种伪随机。
下面介绍srand函数。srand(unsigned seed)通过参数seed改变系统提供的种子值，从而可以使得每次调用rand函数生成的伪随机数序列不同，从而实现真正意义上的“随机”。通常利用系统时间来改变系统的种子值，即srand(time(NULL))，可以为rand函数提供不同的种子值，进而产生不同的随机数序列。注意，time函数包含在time.h头文件内，故声明时需要加入此文件。
具体使用时可以配合RAND_MAX得到想要范围内的随机值。
若要产生信源0，1，则使得生成的随机数除以2取余数即可。

// 信源生成函数
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>
#define data_len XXX
void main()
{srand((unsigned int)time(NULL));int data_beforecoding[data_len];for (int i = 0; i < data_len; i++){data_beforecoding[i] = rand() % 2;}
}

3. 高斯噪声的加入

生成高斯噪声的方式我们采用Box-Muller方法，其描述如下：
如果在(0,1]\left( {0,1} \right](0,1]值域内有两个均匀分布随机变量XXX和YYY，则可以由以下运算得到正态分布随机变量ZZZ
Z=Rcos⁡θZ=R\cos\thetaZ=Rcosθ
(或Z=Rsin⁡θZ=R\sin\thetaZ=Rsinθ) 其中，R=−2ln⁡XR=\sqrt{-2\ln X}R=−2lnX，θ=2πY\theta {\rm{ = 2}}\pi Yθ=2πY。

此时，正态随机变量Z∼N(0,1)Z \sim N\left( {0,1} \right)Z∼N(0,1)。

下面我们来证明一下：
假定随机变量X∼N(0,1)X \sim N\left( {0,1} \right)X∼N(0,1)，Y∼N(0,1)Y \sim N\left( {0,1} \right)Y∼N(0,1)，且相互独立，令p(X)p\left( X \right)p(X)和p(Y)p\left( Y \right)p(Y)分别为其概率密度函数，有：

p(X)=12πe−X22p\left( X \right) = \frac{1}{{\sqrt {2\pi } }}{e^{ - \frac{{{X^2}}}{2}}}p(X)=2π1e−2X2
p(Y)=12πe−Y22p\left( Y \right) = \frac{1}{{\sqrt {2\pi } }}{e^{ - \frac{{{Y^2}}}{2}}}p(Y)=2π1e−2Y2

联合概率密度由于独立，有

p(X,Y)=12πe−X2+Y22p\left( {X,Y} \right) = \frac{1}{{2\pi }}{e^{ - \frac{{{X^2} + {Y^2}}}{2}}}p(X,Y)=2π1e−2X2+Y2

对随机变量XXX和YYY作极坐标变换如下：

X=Rcos⁡θX=R\cos\thetaX=Rcosθ
Y=Rsin⁡θY=R\sin\thetaY=Rsinθ

有：
∫−∞∞∫−∞∞12πe−X2+Y22dXdY=∫02π∫0∞12πe−R22RdθdR=1\int_{ - \infty }^\infty {\int_{ - \infty }^\infty {\frac{1}{{2\pi }}{e^{ - \frac{{{X^2} + {Y^2}}}{2}}}dXdY} } = \int_0^{2\pi } {\int_0^\infty {\frac{1}{{2\pi }}{e^{ - \frac{{{R^2}}}{2}}}Rd\theta dR} } = 1∫−∞∞∫−∞∞2π1e−2X2+Y2dXdY=∫02π∫0∞2π1e−2R2RdθdR=1

由此可得随机变量RRR和θ\thetaθ的分布函数：

PR(R≤r)=∫02π∫0r12πe−R22RdθdR=1−e−r22{P_R}\left( {R \le r} \right) = \int_0^{2\pi } {\int_0^r {\frac{1}{{2\pi }}{e^{ - \frac{{{R^2}}}{2}}}Rd\theta dR} } = 1 - {e^{ - \frac{{{r^2}}}{2}}}PR(R≤r)=∫02π∫0r2π1e−2R2RdθdR=1−e−2r2

Pθ(θ≤ϕ)=∫0ϕ∫0∞12πe−R22RdθdR=ϕ2π{P_\theta }\left( {\theta \le \phi } \right) = \int_0^\phi {\int_0^\infty {\frac{1}{{2\pi }}{e^{ - \frac{{{R^2}}}{2}}}Rd\theta dR} } = \frac{\phi }{{2\pi }}Pθ(θ≤ϕ)=∫0ϕ∫0∞2π1e−2R2RdθdR=2πϕ

此时θ\thetaθ服从[0,2π]\left[ {0,2\pi } \right][0,2π]上的均匀分布，下面对RRR进行变换，使其符合均匀分布，构造

FR(r)=1−e−r22{F_R}\left( r \right) = 1 - {e^{ - \frac{{{r^2}}}{2}}}FR(r)=1−e−2r2

得到反函数：

R=FR−1(Z)=−2ln⁡(1−Z)R = F_R^{ - 1}\left( Z \right) = \sqrt{-2\ln \left( {1 - Z} \right)}R=FR−1(Z)=−2ln(1−Z)

当随机变量Z∼U[0,1]Z \sim U\left[ {0,1} \right]Z∼U[0,1]，RRR的分布函数为FR{F_R}FR，同样的，1−Z∼U[0,1]1-Z \sim U\left[ {0,1} \right]1−Z∼U[0,1]，因此可以选择两个变量U1U_1U1和U2U_2U2，使得

θ=2πU1\theta = 2\pi {U_1}θ=2πU1
1−Z=U21-Z =U_21−Z=U2

即

R=−2ln⁡U2R= \sqrt{-2\ln U_2 }R=−2lnU2

由此带入X=Rcos⁡θX=R\cos\thetaX=Rcosθ，Y=Rsin⁡θY=R\sin\thetaY=Rsinθ，即可证明随机变量X∼N(0,1)X \sim N\left( {0,1} \right)X∼N(0,1)，Y∼N(0,1)Y \sim N\left( {0,1} \right)Y∼N(0,1)。由此可知假设条件正确。

证明完毕后，我们可以根据此定理来构造高斯随机变量：先独立构造分布为U[0,1]U\left[ {0,1} \right]U[0,1]上的变量XXX和YYY；随后进行运算：

Z=−2ln⁡Xsin⁡(2πY)Z = \sqrt { - 2\ln X} \sin \left( {2\pi Y} \right)Z=−2lnXsin(2πY)

则有Z∼N(0,1)Z \sim N\left( {0,1} \right)Z∼N(0,1)。

由此可得构造高斯随机变量的函数：

// Gauss型变量生成函数
double Gaussianrand()
{double U, V, Z;U = (double)(rand() + 1.0) / (double)(RAND_MAX + 1.0);V = (double)(rand() + 1.0) / (double)(RAND_MAX + 1.0);Z = sqrt(-2.0 * log(U)) * sin(2.0 * pi * V);return Z;
}

其中加一（+1.0）是为了防止产生0使得对数运算产生错误；使用墙砖是为了防止整型除以整形得到整形（如3/2=1，3.0/2.0=1.5），具体参见C语言数据类型运算。

4. 信噪比计算

首先应确定，高斯白噪声的功率谱密度，噪声功率与方差的关系为：
P(f)=N02=σ2P\left( f \right) = \frac{{{N_0}}}{2} = {\sigma ^2}P(f)=2N0=σ2
而高斯变量有：
X∼N(0,1)⇔μX+σ∼N(μ,σ2)X \sim N\left( {0,1} \right) \Leftrightarrow \mu X + \sigma \sim N\left( {\mu ,{\sigma ^2}} \right)X∼N(0,1)⇔μX+σ∼N(μ,σ2)
故设信噪比为EbN0_dB（以分贝为单位），首先确定Eb = 1，转化EbN0的单位，有EbN0 = 10 ^ (EbN0_dB / 10)。之后得到N0：
N0=Eb/EbN0=Eb/10EbN0(dB)10{N_0} = {E_b}/{E_b}{N_0} = {E_b}/{10^{\frac{{{E_b}{N_0}\left( {{\rm{dB}}} \right)}}{{10}}}}N0=Eb/EbN0=Eb/1010EbN0(dB)
之后有：
σ=N0/2\sigma {\rm{ = }}\sqrt {{N_0}/2}σ=N0/2
由此可构造符合信噪比要求的高斯白噪声。

三、BPSK仿真结果和理论值对比

显然，判决方式我们选择相干解调。理论误码率为：
Pe=0.5erfc(r){P_e} = 0.5{\mathop{\rm erfc}\nolimits} \left( {\sqrt r } \right)Pe=0.5erfc(r)
其推导参见通信原理内数字通信系统抗噪声性能有关章节。
其中rrr表示码元信噪比E/N0E/{N_0}E/N0（单位不是dB），其与比特信噪比Eb/N0E_b/{N_0}Eb/N0的关系如下：
E/N0=r/kE/{N_0} = r/kE/N0=r/k
其中kkk表示每个码元所包含的比特数，由于调制是BPSK，进制数M=1M=1M=1，故k=1k=1k=1。
经过仿真2000000次，得到的仿真结果和理论值对比如下：

比特信噪比	仿真结果	理论结果
0	0.078697	0.078693
1	0.056417	0.056264
2	0.037619	0.037483
3	0.023006	0.022876
4	0.012448	0.012549
5	0.005925	0.006007
6	0.002422	0.023953
7	0.000755	0.007775
8	0.000196	0.000191
9	0.000039	0.000038
10	0.000000	0.000007

对比可知仿真结果差距不大，当仿真次数增大，按中心极限定理差距将逐渐缩小。

五、BPSK仿真程序代码

代码封装部分省略，这里给出的代码会存在堆栈溢出问题，但解决办法可以通过设置全局变量和函数流程化封装解决。

//BPSK调制 —— C code
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>
#define data_len 1000000   //比特长
#define pi 3.1415926535
#define EbN0_MAX 11     //信噪比 0 - 10double Gaussianrand();
double Erdproduce(double EbN0_dB);int main()
{srand((unsigned int)time(NULL));double EbN0_dB, Erd;for (int i = 0; i < EbN0_MAX; i++){EbN0_dB = (double)i;Erd = Erdproduce(EbN0_dB);printf("%lf \n", Erd);}return 0;
}double Erdproduce(double EbN0_dB)
{int data_beforecoding[data_len];       //信源for (int i = 0; i < data_len; i++){data_beforecoding[i] = rand() % 2;}int data_coding[data_len];             //调制for (int i = 0; i < data_len; i++){if (data_beforecoding[i] == 0)data_coding[i] = 1;elsedata_coding[i] = -1;}double Gaussnoice[data_len];for (int i = 0; i < data_len; i++){Gaussnoice[i] = Gaussianrand();}double Eb = 1;double EbN0;EbN0 = pow(10, EbN0_dB / 10.0);double N0 = Eb / EbN0;double sigma = sqrt(0.5 * N0);for (int i = 0; i < data_len; i++){Gaussnoice[i] = Gaussnoice[i] * sigma;}//加入高斯噪声double re_symbol[data_len];for (int i = 0; i < data_len; i++){re_symbol[i] = data_coding[i] + Gaussnoice[i];}//判决int re_data[data_len];for (int i = 0; i < data_len; i++){if (re_symbol[i] < 0)re_data[i] = 1;elsere_data[i] = 0;}//误比特统计double Erd = 0.0;int Nrd = 0;for (int i = 0; i < data_len; i++){if (re_data[i] != data_beforecoding[i])Nrd = Nrd + 1;}Erd = (double)Nrd / (double)data_len;return Erd;
}
//高斯白噪声子函数
double Gaussianrand()
{double U, V, Z;U = (double)(rand() + 1.0) / (double)(RAND_MAX + 1.0) ;V = (double)(rand() + 1.0) / (double)(RAND_MAX + 1.0) ;Z = sqrt(-2.0 * log(U)) * sin(2.0 * pi * V);return Z;
}

六、QPSK仿真程序

QPSK是经过封装后的程序，避免了溢出问题。

#include<stdio.h>
#include<stdlib.h>
#include<time.h>
#include<math.h>
#define pi 3.1415926535
#define data_len 1000000
#define EbN0_dB_MAX 11
#define k 2                 //符号比特数
#define M 2                 //进制数const int symbol_len = data_len / k;
int EbN0_dB;
const int Eb = 1;
int data_bfcoding[data_len], data_decoding[data_len];
int symbol_bfcoding[data_len];
int symbol_I[symbol_len], symbol_Q[symbol_len];
double symbol_I_AWGN[symbol_len], symbol_Q_AWGN[symbol_len];
int symbol_I_decoding[symbol_len], symbol_Q_decoding[symbol_len];
void Data_generate();
void Symbol_generate();
void Modulation();
double Gaussnoise();
void AWGN(int EbN0_dB);
void Demodulation();
double Erdsimulation();int main()
{int j = 0;double Erd; srand((unsigned int(time(NULL))));for (int j = 0; j < EbN0_dB_MAX; j++){EbN0_dB = j;Data_generate();Symbol_generate();Modulation();AWGN(EbN0_dB);Demodulation();Erd = Erdsimulation();printf("%f \n", Erd);}return 0;}void Data_generate()
{for (int i = 0; i < data_len; i++){data_bfcoding[i] = rand() % 2;}
}void Symbol_generate()
{for (int i = 0; i < data_len; i++){if (i % 2 == 0)symbol_I[i / 2] = data_bfcoding[i];elsesymbol_Q[i / 2] = data_bfcoding[i];}
}void Modulation()
{for (int i = 0; i < symbol_len; i++){symbol_I[i] = 1 - 2 * symbol_I[i];symbol_Q[i] = 1 - 2 * symbol_Q[i];}
}double Gaussnoise()
{double U, V, Z;U = (double)(rand() + 1.0) / (double)(RAND_MAX + 1.0);V = (double)(rand() + 1.0) / (double)(RAND_MAX + 1.0);Z = sqrt(-2 * log(U)) * sin(2 * pi * V);return Z;
}void AWGN(int EbN0_dB)
{double EbN0 = pow(10, (double)EbN0_dB / 10);double N0 = Eb / EbN0;double sigma = sqrt(0.5 * N0);for (int i = 0; i < symbol_len; i++){symbol_I_AWGN[i] = (double)symbol_I[i] + Gaussnoise() * sigma;symbol_Q_AWGN[i] = (double)symbol_Q[i] + Gaussnoise() * sigma;}
}void Demodulation()
{for (int i = 0; i < symbol_len; i++){if (symbol_I_AWGN[i] < 0.0){symbol_I_decoding[i] = 1;}else{symbol_I_decoding[i] = 0;}if (symbol_Q_AWGN[i] < 0.0){symbol_Q_decoding[i] = 1;}else{symbol_Q_decoding[i] = 0;}}for (int i = 0; i < data_len; i++){if (i % 2 == 0)data_decoding[i] = symbol_I_decoding[i / 2];elsedata_decoding[i] = symbol_Q_decoding[i / 2];}}double Erdsimulation()
{int Nrd = 0;double Erd_temp;for (int i = 0; i < data_len; i++){if (data_decoding[i] != data_bfcoding[i])Nrd++;}Erd_temp = (double)Nrd / (double)data_len;return Erd_temp;
}

七、参考资料

BPSK调制与解调-MATLAB基带仿真
C语言随机数：rand()和srand(time(NULL))的使用
Box-Muller变换原理详解