对于调制信号，我们可以通过其功率谱（PSD）来研究其频谱特性。因为无线电频谱是有限且十分珍贵的资源，一个调制信号的 PSD 能够反应其频谱占用情况，是研究信号特性的重要工具。

# 随机过程

随机过程是通信系统中的一个基础概念，也是通信系统中的一个基础概念。

一个随机过程$X(t)$ 是一组依赖于变量$t$ 的随机变量。定义随机过程$X(t)$ 的均值$m_X(t)$ 和自相关函数$R_X(t_1,t_2)$

$$m_X(t) = \mathbb{E}[X(t)]\\[2ex] R_X(t_1,t_2) = \mathbb{E}[X(t_1)X^*(t_2)]$$

两个随机过程$X(t)$ 和$Y(t)$ 的互相关函数

$$R_{XY}(t_1,t_2)= \mathbb{E}[X(t_1)Y^*(t_2)]$$

如果$R_X(t_1,t_2)=R_X^*(t_1,t_2)$ 那么这个过程是 Hermitian 对称的，互相关函数定义类似。

如果一个随机过程$X(t)$ 的均值是一个常数，$R_X(t_1,t_2)=R_X(\tau)$，$\tau=t_1-t_2$，这个过程就是广义平稳过程（WSS）。WSS 随机过程的功率谱密度是描述了其功率关于频率分布的函数$S_X(f)$。Wiener-Khinchin 定理说明了，WSS 过程的功率谱是自相关函数$R_X(\tau)$ 的 Fourier 变换：

$$\mathcal{S}_X(f)=\mathcal{F}[R_X(\tau)]$$

如果随机过程$X(t)$ 的均值和自相关函数都是以$T_0$ 周期的周期函数，那么这个随机过程是一个循环平稳（Cyclostationary）过程。循环平稳过程满足：

$$m_X(t+T_0)=m_X(t)\\[2ex] R_X(t_1+T_0,t_2+T_0)=R_X(t_1,t_2)$$

对于循环平稳过程，其平均自相关函数是其在一个周期内的自相关函数的平均值：

$$\overline{R_X(t)}\frac{1}{T_0}\int_{0}^{T_0}R_X(t+\tau,t){\rm d}t$$

循环平稳随机过程的平均功率谱密度为平均自相关函数的 Fourier 变换：

$$\overline{S_X(f)}=\mathcal{F}[\overline{R_X(t)}]$$

# 一般信号的功率谱

我们写一个一般情况下的低通信号：

$$\boldsymbol{v}_1(t)=\sum_{n=-\infty}^{\infty}s(t-nT;\boldsymbol{I_n})$$

那么这个信号的数学期望是

$$\mathbb{E}[\boldsymbol{v}_k(t)]=\sum_{n=-\infty}^{\infty}\mathbb{E}[\boldsymbol{I_n}]g(t-nT)=\mu_{\boldsymbol{I}}\sum_{n=-\infty}^{\infty}g(t-nT)=\mathbb{E}[\boldsymbol{v}_1(t+T)]$$

这个信号的自相关函数

$$R_{v1}(t_1,t_2)=\mathbb{E}[\boldsymbol{v}_1(t_1)\boldsymbol{v}_1^*(t_2)]=\sum_{n=-\infty}^{\infty}\sum_{m=-\infty}^{\infty}\mathbb{E}[\boldsymbol{I}_n\boldsymbol{I}_m^*]g(t_1-nT)g^*(t_2-mT)=R_{v1}(t_1+T,t_2+T)$$

上面的数学期望和自相关函数说明了这个信号是一个循环平稳过程。因此，我们可以通过平均自相关函数的 Fourier 变换来计算其功率谱密度：

$$\begin{aligned} &\overline{R}_{v1}(t)\\ =&\frac{1}{T}\int_0^T\sum_{n=-\infty}^{\infty}\sum_{m=-\infty}^{\infty}R_{\boldsymbol{I}}(n-m)g(t+\tau-nT)g^*(t+\tau-mT){\rm d}\tau\\ \stackrel{k=n-m}{=}&\frac{1}{T}\sum_{k=-\infty}^{\infty}R_{\boldsymbol{I}}(k)\sum_{m=-\infty}^{\infty}\int_0^Tg(t+\tau-kT-mT)g^*(t-mT){\rm d}t\\ \stackrel{u=t-mT}{=}&\frac{1}{T}\sum_{k=-\infty}^{\infty}R_{\boldsymbol{I}}(k)\sum_{m=-\infty}^{\infty}\int_{-mT}^{-(m-1)T}g(u+\tau-kT)g^*(u){\rm d}u\\ =&\frac{1}{T}\sum_{k=-\infty}^{\infty}R_{\boldsymbol{I}}(k)\int_{-\infty}^{\infty}g(u+\tau-kT)g^*(u){\rm d}u \end{aligned}$$

令

$$g_m(\tau)=R_{\boldsymbol{I}}(m)\int_{-\infty}^{\infty} g(u+\tau)g^*(u){\rm d}u$$

带入得到

$$\begin{aligned} &\overline{R}_{v1}(t)\\ =&\frac{1}{T}\sum_{k=-\infty}^{\infty}g_m(\tau-kT) \end{aligned}$$

这部分式子依旧是十分复杂的，难以进行计算。然而，我们观察$g_m(\tau)$，应该能够发现这个函数的形式类似卷积。考虑到卷积的性质，这个函数在频域上的表达式应该会比较清晰。

我们将$g_m(\tau)$ 移到频域上，得到：

$$\begin{aligned} &G_m(f)\\ =&\int_{-\infty}^{\infty}g_m(\tau)e^{-j2\pi f\tau}{\rm d}\tau\\ =&\int_{-\infty}^{\infty}(R_{\boldsymbol{I}}(m)\int_{-\infty}^{\infty} g(u+\tau)g^*(u){\rm d}u)e^{-j2\pi f\tau}{\rm d}\tau\\ =&R_{\boldsymbol{I}}(m)\int_{-\infty}^{\infty}(\int_{-\infty}^{\infty}g(u+\tau)e^{-j2\pi f\tau}{\rm d}\tau)g^*(u){\rm d}u\\ \stackrel{v=u+t}{=}&R_{\boldsymbol{I}}(m)\int_{-\infty}^{\infty}(\int_{-\infty}^{\infty}g(v)e^{-j2\pi fv}{\rm d}v)g^{*}(u){\rm d}u\\ =&R_{\boldsymbol{I}}(m)\int^{\infty}_{-\infty}g(v)e^{-j2\pi fv}{\rm d}v\int^{\infty}_{-\infty}g^{*}(u)e^{-j2\pi fu}{\rm d}u\\ =&R_{\boldsymbol{I}}(m)|G(f)|^2 \end{aligned}$$

这样子，对上面自相关函数进行 Fourier 变换，并且带入这里的$G_m(f)$，我们就得到了这个信号的功率谱密度函数：

$$\begin{aligned} &\overline{S}_{v1}(f)\\ =&\mathcal{F}\{\overline{R}_{v1}(t)\}\\ =&\frac{1}{T}\sum_{k=-\infty}^{\infty}\mathcal{F}\{g_k(t-kT)\}\\ =&\frac{1}{T}\sum_{k=-\infty}^{\infty}R_{\boldsymbol{I}}(m)|G(f)|^2e^{-2\pi kfT}\\ =&\frac{1}{T}S_{\boldsymbol{I}}(f)|G(f)|^2 \end{aligned}$$

其中

$$S_{I}(f)=\sum_{k=-\infty}^{\infty}R_I(k)e^{-j2\pi kfT}$$

这样子，我们就得到了数字通信的功率谱的一般表达式。这个表达式说明了两个因素决定了功率谱的形状，一个是信息序列的自相关函数，另一个是 pulse shaping function 的频谱。

# 互不相关的信息序列下的功率谱

考虑输入一个互不相关的信息序列，其自相关函数为

$$R_{I}(k)=\begin{cases} \mu_I^2+\sigma_I^2&,k=0\\ \mu_I^2&,k\neq 0 \end{cases}$$

那么

$$S_I(f)=\sigma_I^2+\mu_I^2\sum_{l=-\infty}^{\infty}e^{-j2\pi fkT}$$

根据 Poisson 求和公式：

$$\sum_{k=-\infty}^{\infty}e^{-j2\pi fkT}=\frac{1}{T}\sum_{-\infty}^{\infty}\delta(f-\frac{k}{T})$$

带入$S_I(f)$

$$S_I(f)=\sigma_I^2+\frac{\mu_I^2}{T}\sum_{k=-\infty}^{\infty}\delta(f-\frac{k}{T})$$

带入到$\overline{S}_{v1}(f)$ 中，我们得到：

$$\overline{S}_{v1}(f)=\underbrace{\frac{\sigma_I^2}{T}|G(F)|^2}_{continuous}+\underbrace{\frac{\mu_I^2}{T}\sum_{k=-\infty}^{\infty}\delta(f-\frac{k}{T})|G(f)|^2}_{discrete}$$

也就是这样的一个信号的功率谱分为了连续和离散的两个部分，如果信息序列的数学期望为 0，那么这个信号的功率谱就没有离散的部分，功率谱的形状完全由脉冲整形函数决定。在实际的使用中，我们一般是希望有尽可能少的冲激函数部分，因此我们都会尽可能希望信息序列的数学期望为 0。

# 使用方波的情况

如果我们使用方波作为$g(t)$，即

$$g(t)=A[u_{-1}-u_{-1}(t-T)]$$

其 Fourier 变换为

$$G(F)=AT{\rm sinc}(fT)e^{-j\pi fT}$$

那么其功率谱为

$$\begin{aligned} &\overline{S_v1}(f)\\ =&\frac{\sigma_I^2}{T}|G(f)|^2+\frac{\mu_I^2}{T}\sum_{k=-\infty}^{\infty}\delta(f-\frac{k}{T})|G(f)|^2\\ =&\sigma^2_IA^2T{\rm sinc}^2(fT)+\mu_I^2A^2\delta(f) \end{aligned}$$

我们知道${\rm sinc}(x)$ 函数有一个特殊的性质，对于$x=\pm 1,\pm 2,...$，${\rm sinc}(x)$ 的值为 0，这样就使得在$\mu_I\neq 0$ 的情况下，信号的功率谱也就只会在$f=0$ 处有一个高度为$\mu_I^2A^2$ 的冲激函数。

# 使用升余弦函数的情况

如果我们使用升余弦函数作为$g(t)$，即

$$g(t)=\frac{A}{2}[1+\cos(\frac{2\pi}{T}(t-\frac{T}{2}))](u_{-1}(t)-u_{-1}(t-T))$$

其 Fourier 变换为

$$G(f)=\frac{AT}{2}{\rm sinc}(fT)\frac{1}{1-f^2T^2}e^{-j\pi fT}$$

带入其功率谱，得到

$$\overline{S_v1}(f) =\frac{\sigma_I^2A^2T{\rm sinc}^2(fT)}{4(1-f^2T^2)^2}+\frac{\mu_I^2A^2}{4}\delta(f)+\frac{\mu_I^2A^2}{16}\delta(f-\frac{1}{T})+\frac{\mu_I^2A^2}{16}\delta(f+\frac{1}{T})$$

这里注意其中一步化简需要用到一个常用的极限：

$$\lim_{x\to \pm 1}\frac{\text{\rm sinc}^2(x)}{(1-x^2)^2}=\frac{1}{4}$$

如果从$\overline{S}_{v1}(f)$ 来看的话，使用 RC 函数的功率谱关于$f$ 的衰减为$1/f^{6}$，而使用了方波的情况为$1/f^2$，这也就意味着使用 RC 函数的情况下，信号的能量更加集中。

我们绘制两种情况下的功率谱图：

可以看到 RC 的旁瓣更大，但是衰减速度更快，能量更加集中在中心频率上。一般来说，脉冲整形函数的波形越平滑（能够连续微分的阶数更高），其功率谱衰减更快。

当然，尽管 RC 的性能是优于方波的，但是其实现成本更高，这也需要进行取舍。

# 相互关联的信息序列的功率谱

考虑一个特殊的情况，信息序列为$I_n=b_n+b_{n+1}$，其中$b_n$ 为互不相关，均值为 0，方差为 1 的信息序列。此时$I_n$ 的自相关函数：

$$R_I(k)=\begin{cases} 2&,k=0\\ 1&,k=1\\ 0&,otherwise \end{cases}$$

此时信息序列的功率谱为：

$$\begin{aligned} S_I(f)=2+B^{j2\pi fT}+e^{-j2\pi fT}=4\cos^2(\pi fT) \end{aligned}$$

这样频谱的宽度更窄，同时衰减速度也变快。

# CPFSK 和 CPM 的功率谱

CPM 信号的功率谱推导过程和前面的稍有区别，我们考虑独立同分布的信息序列$\boldsymbol{I}$，CPM 信号表示为

$$\boldsymbol{v}_1(t)=e^{j\phi(t;\boldsymbol{I})}$$

其中

$$\phi(t;\boldsymbol{I})=2\pi h\sum_{k=-\infty}^{\infty}I_kq(t-kT)$$

我们计算其自相关函数

$$\begin{aligned} &R_{v1}(t_1,t_2)\\ =&\mathbb{E}[\boldsymbol{v}_1(t_1)\boldsymbol{v}^*_1(t_2)]\\ =&\mathbb{E}[e^{j\phi(t_1;\boldsymbol{I})}e^{-j\phi(t_2;\boldsymbol{I})}]\\ =&\mathbb{E}[\exp(j2\pi h\sum_{k=-\infty}^{\infty}I_k[q(t_1-kT)-q(t_2-kT)])]\\ =&\mathbb{E}[\prod_{k=-\infty}^{\infty}\exp(j2\pi hI_k[q(t_1-kT)-q(t_2-kT)])]\\ =&\prod_{k=-\infty}^{\infty}\mathbb{E}[\exp(j2\pi hI_k[q(t_1-kT)-q(t_2-kT)])]\\ =&\prod_{k=-\infty}^{\infty}[\sum_{n\in\mathcal{S}}P_n\exp(j2\pi hn[q(t_1-kT)-q(t_2-kT)])] \end{aligned}$$

这里记所有可能取值的码字$I_k=n\in\mathcal{S}$，其中$I_k=n$ 的概率为$P_n$。上述步骤中一些转换需要利用到独立同分布的性质。

接下来，我们计算周期内的平均自相关函数：

$$\begin{align} &\overline{R}_{v1}(\tau)\\ =&\frac{1}{T}\int_0^T R_{v1}(t+\tau,t){\rm d}t\\ =&\frac{1}{T}\int_0^T\prod_{k=-\infty}^{\infty}[\sum_{n\in\mathcal{S}}P_n\exp(j2\pi hn[q(t+\tau-kT)-q(t-kT)])]{\rm d}t\\ \end{align}$$

这个很不好，包含了积分连乘和累加，同时连乘还是无限多项的。但是我们注意到，如果指数中的$j2\pi hn[q(t+\tau-kT)-q(t-kT)])=0$，那么这一部分累乘中间的部分就是$\sum_{n\in \mathcal{S}}P_ne^0=1$，这部分就不用参加累乘计算了。

我们那么就需要来考量一下这个部分的特性。

前面一节提到了 phase shaping function $q(t)$ 的性质，这一部分在$t<0$ 时值为 0，在$t>LT$ 时值为$1/2$。

我们设置一个整数$m$，满足$mT\le\tau=\xi + mT<(m+1)T$，$\xi\in[0,1)$。

这样子对于$q(t-kT)$ 和$q(t+\tau-kT)$，$t\in[0,T)$，在$k<\min\{1-L,m+1-L\}$ 时，两个取值都是 0，在$k>\max\{0,m+1\}$ 时，两个取值都是$1/2$，这两段区间上两个函数相减，指数均为 0，可以不用讨论。

那么平均自相关函数就变成了：

$$\begin{align} &\overline{R}_{v1}(\tau)\\ =&\frac{1}{T}\int_{0}^{T}\prod_{k=\min\{m+1-L,1-L\}}^{\max\{m+1,0\}}[\sum_{n\in\mathcal{S}}P_ne^{j2\pi hn[q(t+\tau-kT)-q(t-kT)]}]{\rm d}t\\ =&\frac{1}{T}\int_{0}^{T}\prod_{k=1-L}^{m+1}[\sum_{n\in\mathcal{S}}P_ne^{j2\pi hn[q(t+\tau-kT)-q(t-kT)]}{\rm d}t \end{align}$$

接下来把平均自相关函数进行 Fourier 变换，得到功率谱密度。这里补充一点，这里的平均自相关函数满足 Hermitian 对称性，所以

$$\overline{R}_{v1}(-\tau)=\overline{R}_{v1}^*(\tau)$$

进行 Fourier 变换：

$$\begin{aligned} &\overline{S}_{v1}(f)\\ =&\int_{-\infty}^{\infty}\overline{R}_{v1}(\tau)e^{-j2\pi f\tau}{\rm d}\tau\\ =&\int_{-\infty}^{0}\overline{R}_{v1}(\tau)e^{-j2\pi f\tau}{\rm d}\tau+\int_{0}^{\infty}\overline{R}_{v1}(\tau)e^{-j2\pi f\tau}{\rm d}\tau\\ =&\int_{0}^{\infty}\overline{R}_{v1}(-\tau)e^{j2\pi f\tau}{\rm d}\tau+\int_{0}^{\infty}\overline{R}_{v1}(\tau)e^{-j2\pi f\tau}{\rm d}\tau\\ =&\int_{0}^{\infty}(\overline{R}_{v1}(\tau)e^{-j2\pi f\tau})^*{\rm d}\tau+\int_{0}^{\infty}\overline{R}_{v1}(\tau)e^{-j2\pi f\tau}{\rm d}\tau\\ =&2{\rm Re}\{\int_0^{\infty}\overline{R}_{v1}(\tau)e^{-j2\pi f\tau}{\rm d}\tau\} \end{aligned}$$

这里面的$\int_0^{\infty}\overline{R}_{v1}(\tau)e^{-j2\pi f\tau}{\rm d}\tau$ 展开，这一部分中，我们注意到，如果$m\ge L$，对于$q(t-kT)$ 和$q(t+\tau-kT)$，$t\in[0,T)$，$m+1-L>0$，此时两个 phase shaping function 变化段没有重叠！

$$\begin{aligned} &\int_0^{\infty}\overline{R}_{v1}(\tau)e^{-j2\pi f\tau}{\rm d}\tau\\ =&\int_{0}^{LT}\overline{R}_{v1}(\tau)e^{-j2\pi f\tau}{\rm d}\tau+\int_{LT}^{\infty}\overline{R}_{v1}(\tau)e^{-j2\pi f\tau}{\rm d}\tau\\ =&\int_{0}^{LT}\overline{R}_{v1}(\tau)e^{-j2\pi f\tau}{\rm d}\tau+\sum_{m=L}^{\infty}\int_{mT}^{(m+1)T}\overline{R}_{v1}(\tau)e^{-j2\pi f\tau}{\rm d}\tau \end{aligned}$$

这时候前面一项范围有限，但是后面一项还有无穷级数，不过上面说了，对于$m\ge L$ 的情况，两个 phase shaping function 变化段没有重叠，因此：

$$\begin{aligned} &\overline{R}_{v1}(\tau)\\ =&\frac{1}{T}\int_{0}^{T}\prod_{k=1-L}^{m+1}[\sum_{n\in S}P_ne^{j2\pi hn[q(t+\tau-kT)-q(t-kT)]}]{\rm d}t\\ =&\frac{1}{T}([\prod_{k=1-L}^{0}P_ne^{j2\pi hn[q(t+\tau-kT)-q(t-kT)]}]\\ &\times[\prod_{k=1}^{m-L}P_ne^{j2\pi hn[q(t+\tau-kT)-q(t-kT)]}]\\ &\times[\prod_{k=M-L+1}^{m+1}P_ne^{j2\pi hn[q(t+\tau-kT)-q(t-kT)]}])\\ \end{aligned}$$

这里将连乘中的项分为了三个部分，这三个部分分别对应了$q(t-kT)$ 和$q(t+\tau-kT)$ 各自的变化段和中间的不变段，因此我们得到：

$$\begin{aligned} &\overline{R}_{v1}(\tau)\;(m\ge L)\\ =&\frac{1}{T}([\prod_{k=1-L}^{0}P_ne^{j2\pi hn[\frac{1}{2}-q(t-kT)]}]\\ &\times[\prod_{k=1}^{m-L}P_ne^{j2\pi hn[\frac{1}{2}-0]}]\\ &\times[\prod_{k=M-L+1}^{m+1}P_ne^{j2\pi hn[q(t+\tau-kT)-0]}])\\ \end{aligned}$$

令

$$\phi_I(h)=\sum_{n\in \mathcal S}P_ne^{j\pi hn}$$

称为特征函数，则有：

$$\begin{aligned} &\overline{R}_{v1}(\tau)\;(m\ge L)\\ =&\frac{1}{T}\int_{0}^{T}(\prod_{k=1-L}^{0}[\sum_{n\in \mathcal S}P_ne^{j2\pi hn[\frac{1}{2}-q(t-kT)]}][\phi_I(h)]^{m-L})[\prod_{k'=1-L}^{1}[\sum_{n\in \mathcal S}P_ne^{j2\pi hn[q(t+\tau-k'T-mT)]}]]{\rm d}t\\ =&[\phi_I(h)]^{m-L}\lambda(\tau-mT) \end{aligned}$$

这里$\lambda(\xi)$ 定义为：

$$\lambda(\xi)=\frac{1}{T}\int_{0}^{T}(\prod_{k=1-L}^{0}[\sum_{n\in\mathcal S}P_ne^{j2\pi hn[\frac{1}{2}-q(t-kT)]}]\prod_{k'=1-L}^{1}[\sum_{h\in\mathcal S}P_ne^{j2\pi hnq(t+\xi-k'T)}]){\rm d}t$$

继续带入前面的$\int_0^{\infty}\overline{R}_{v1}(\tau)e^{-j2\pi f\tau}{\rm d}\tau$，得到

$$\begin{aligned} &\int_0^{\infty}\overline{R}_{v1}(\tau)e^{-j2\pi f\tau}{\rm d}\tau\\ =&\sum_{m=L}^{\infty}\int_{mT}^{(m+1)T}[\phi_I(h)]^{m-L}\lambda(\tau-mT)e^{-j2\pi f\tau}{\rm d}\tau\\ =&\sum_{m=L}^{\infty}\int_{0}^{T}[\phi_I(h)]^{m-L}\lambda(\xi)e^{-j2\pi f(\xi+mT)}{\rm d}\xi,\;(\xi=\tau-mT)\\ =&({\color{blue}\sum_{m=L}^{\infty}\phi_I(h)^{m-L}})（\int_{0}^{T}\lambda(\xi)e^{-j2\pi f\xi}{\rm d}\xi)\\ =&\begin{cases} \frac{e^{-j\pi fLT}}{1-\phi_I(h)e^{-j2\pi fT}}\int_{0}^{T}\lambda(\xi)e^{-j2\pi f\xi}{\rm d}\xi,&|\phi_I(h)|<1\\[2ex] (e^{-j2\pi fLT}\sum_{m'=0}^\infty e^{-j2\pi T(f-v/T)m'})(\int_0^\infty\lambda(\xi)e^{-j2\pi f\xi}{\rm d}\xi),&|\phi_I(h)|=|e^{j2\pi v}|=1 \end{cases} \end{aligned}$$

这边蓝色部分是一个复数的等比数列，我们能够套用求和公式。不过这里有个问题，公比可以是 1，这就不好办了。在这种情况下，我们令

$$\phi_I(h)=e^{j2\pi v}$$

下面来看一下这个等比数列特殊情况下的求和。前面提到了 Poisson 求和公式，我们这里也要用到。现在考虑一个单位阶跃函数：

$$u_{-1}(t)=\begin{cases} 1,&t>0\\ \frac{1}{2},&t=0\\ 0,&t<0 \end{cases}$$

其 Fourier 变换为

$$U_{-1}(f)=\frac{1}{2}(\delta(f)-j\frac{1}{\pi f})$$

我们通过间隔$T_s$ 无限冲激序列对这个函数进行采样，利用冲激函数的性质：

$$\begin{aligned} &U_{-1,\delta}(f)\\ =&\sum_{n=-\infty}^{\infty}u_{-1}(nT_s)e^{-j2\pi nT_sf}\\ =&\frac{1}{2}+\sum_{n=1}^{\infty}e^{-j2\pi nT_sf}\\ =&-\frac{1}{2}+\sum_{n=0}^{\infty}e^{-j2\pi nT_sf} \end{aligned}$$

而套用 Poisson 求和公式，我们得到：

$$\begin{aligned} &U_{-1,\delta}(f)\\ =&\frac{1}{T_s}\sum_{n=-\infty}^{\infty}U_{-1}(f-\frac{n}{T_s})\\ =&\frac{1}{2T_s}\sum_{n=-\infty}^{\infty}(\delta(f-\frac{n}{T_s}-j\frac{1}{\pi(f-\frac{n}{T_s})})) \end{aligned}$$

联立两式，我们得到

$$\sum_{m'=0}^\infty e^{-j2\pi T(f-v/T)m'}=\frac{1}{2}+\frac{1}{2T}\sum_{m'=-\infty}^{\infty}(\delta(f-\frac{v+m'}{T})-j\frac{1}{\pi(f-\frac{v+m'}{T})})$$

这样，我们能够解决前面等比数列求和的问题了：

$$\begin{aligned} &\int_0^{\infty}\overline{R}_{v1}(\tau)e^{-j2\pi f\tau}{\rm d}\tau\\ =&\begin{cases} \frac{e^{-j\pi fLT}}{1-\phi_I(h)e^{-j2\pi fT}}\int_{0}^{T}\lambda(\xi)e^{-j2\pi f\xi}{\rm d}\xi,&|\phi_I(h)|<1\\[2ex] (e^{-j2\pi fLT}\sum_{m'=0}^\infty e^{-j2\pi T(f-v/T)m'})(\int_0^\infty\lambda(\xi)e^{-j2\pi f\xi}{\rm d}\xi),&|\phi_I(h)|=|e^{j2\pi v}|=1 \end{cases}\\ =&\begin{cases} \frac{e^{-j\pi fLT}}{1-\phi_I(h)e^{-j2\pi fT}}\int_{0}^{T}\lambda(\xi)e^{-j2\pi f\xi}{\rm d}\xi,&|\phi_I(h)|<1\\[2ex] e^{-j2\pi fLT}(\frac{1}{2}+\frac{1}{2T}\sum_{m'=-\infty}^{\infty}(\delta(f-\frac{v+m'}{T})-j\frac{1}{\pi(f-\frac{v+m'}{T})}))(\int_0^\infty\lambda(\xi)e^{-j2\pi f\xi}{\rm d}\xi),&|\phi_I(h)|=|e^{j2\pi v}|=1 \end{cases} \end{aligned}$$

现在，对于$|\phi_I(h)|<1$ 的情况，我们可以计算：

$$\begin{aligned} &\overline{S}_{v1}(f) =&2{\rm Re}\{\int_0^{LT}\overline{R}_{v1}(\tau)e^{-j2\pi f\tau}{\rm d}\tau+\frac{1}{1-\phi_I(h)e^{-j2\pi fT}}(\int_0^T\lambda({\xi})e^{-j2\pi f(\xi+LT)}{\rm d}\xi)\} \end{aligned}$$

其中对于$0\le\tau=\xi+mT\le LT$

$$\overline{R}_{v1}(\tau)=\frac{1}{T}\int_{0}^{T}\prod_{k=1-L}^{m+1}[\sum_{n\in\mathcal S}P_ne^{j2\pi hn[q(t+\tau-kT)-q(t-kT)]}]{\rm d}t$$

而在$|\phi_I(h)|=1$ 时，在$f_{m'}=\frac{v+m'}{T}$ 出存在冲激函数。

而如果$|\phi_I(h)|<1$，在符号等概率的情况下：

$$P_n=\frac{1}{M}$$

特征函数

$$\phi_I(h)=\frac{1}{M}\sum_{n=-(M-1)}^{M-1}e^{j\pi hn}=\frac{1}{M}\frac{\sin M\pi h}{\sin \pi h}$$

简化平均自相关函数

$$\overline{R}_{v1}(\tau)=\frac{1}{2T}\int_0^T\prod_{k=1-L}^{\tau/T}\frac{1}{M}\frac{\sin(2\pi hM[q(t+\tau-kT)-q(t-kT)])}{\sin(2\pi h[q(t+\tau-kT)-q(t-kT)])}{\rm d}t$$

相应的功率谱密度

$$S_{v1}(f)=2[\int_0^{LT}\overline{R}_{v1}(\tau)\cos(2\pi f\tau){\rm d}\tau\\+\frac{1-\phi_I(h)\cos(2\pi fT)}{1+\phi_1^2(h)-2\phi_I(h)\cos(2\pi fT)}\int_{LT}^{(L+1)T}\overline{R}_{v1}(\tau)\cos(2\pi f\tau){\rm d}\tau\\-\frac{1-\phi_I(h)\sin(2\pi fT)}{1+\phi_1^2(h)-2\phi_I(h)\sin(2\pi fT)}\int_{LT}^{(L+1)T}\overline{R}_{v1}(\tau)\sin(2\pi f\tau){\rm d}\tau]$$

说实话上面这个公式真套进去跑数值计算也非常费劲…… 两层积分套在一起实在没办法搞……

# CFPSK 的 PSD

如果考虑到使用方波的情况，线性的$q(t)$ 会让整个式子变得比较好看。这里直接拿书上的公式：

$$S_v(f)=T[\frac{1}{M}\sum_{n=1}^{M}A_n^2(f)+\frac{2}{M^2}\sum_{n=1}^{M}\sum_{m=1}^{M}B_{nm}(f)A_n(f)A_m(f)]\\ A_n(f)=\frac{\sin\pi[fT-\frac{1}{2}(2n-1-M)h]}{\pi[fT-\frac{1}{2}(2n-1-M)h]}\\ B_{nm}(f)=\frac{\cos(2\pi fT-\alpha_{nm})-\phi_I(h)\cos\alpha_{nm}}{1+\phi_I^2(h)-2\phi_I\cos(2\pi fT)}\\ \alpha_{nm}=\pi h(m+n-1-M)$$

对于$M=2$ 的情况，我们可以得到：

可以注意到，当$h$ 趋近于 1 时，PSD 中出现了尖峰。在$h$ 小于 1 时，能量主要集中在频谱中心部分，随着$h$ 逐渐趋向于 1，频谱有向外产生一个尖峰的趋势。对于$h$ 大于 1 的情况，PSD 能量主要集中在外侧一般区域中，随着$h$ 逐渐趋向于 1，频谱有向内的趋势。

为了尽可能减少频谱占用，一般我们还是要把调制指数设置在$1$ 以内。当然，在一些特殊的情况下，也许可以通过大于 1 的调制指数来达到一些特殊的效果。

对于$M$ 更高的情况，有类似的：

绘制$M=4$ 和$M=8$，和前面的情况类似。

下面考虑一些特殊情况，我们绘制 MSK 和 OQPSK 的功率谱，能够得到：

可见 MSK 的旁瓣功率下降远快于 OQPSK。

# 附 - 绘图代码

import matplotlib.pyplot as plt

import numpy as np

def sinc(x):

    if x == 0:

        return 1

    return np.sin(np.pi * x) / (np.pi * x)

def s_rec(sigma, A, T, f):

    return pow(sigma, 2) * pow(A, 2) * T * pow(sinc(f * T), 2)

def s_rc(sigma, A, T, f):

    return pow(sigma, 2) * pow(A, 2) * T * pow(sinc(f * T), 2) / (4 * pow((1 - pow(f, 2) * pow(T, 2)),2))

sigma = 1

T = 2

A = 1

x = np.linspace(-10, 10, 10000)

y1 = []

y2 = []

y1_0 = s_rec(sigma, 1, T, 0)

y2_0 = s_rc(sigma, 1, T, 0)

for i in x:

    y1.append(s_rec(sigma, 1, T, i))

    y2.append(s_rc(sigma, 1, T, i))

y1 = np.array(y1)

y2 = np.array(y2)

plt.suptitle("$PSD\ of\ RC\ and\ REC\ (A=1, T=2, \sigma=1)$")

ax = plt.subplot(211)

ax.set_title("$PSD\ of\ RC\ and\ REC$")

plt.plot(x, y1, label="REC")

plt.plot(x, y2, label="RC")

plt.legend()

plt.grid(True)

ax.set_xlabel("$f$")

ax.set_ylabel("$PSD$")

ax = plt.subplot(212)

ax.set_title("$Decay\ RC\ and\ REC\ (\log 10)$")

plt.plot(x, np.log10(y1 / y1_0), label="REC")

plt.plot(x, np.log10(y2 / y2_0), label="RC")

plt.legend()

plt.grid(True)

ax.set_xlabel("$f$")

ax.set_ylabel("$\\log_{10}\\frac{\\overline{S}_{v1}(f)}{\\overline{S}_{v1}(0)}$")

plt.show()

import matplotlib.pyplot as plt

import numpy as np

def phi_h(h, M):

    return np.sin(M * h * np.pi) / np.sin(h * np.pi) / M

def alpha(n, m, M, h):

    return np.pi * (m + n - 1 - M) * h

def A(n, f, T, M, h):

    return (np.sin(np.pi * (f * T - 0.5 * (2 * n - 1 - M) * h))) / (

        np.pi * (f * T - 0.5 * (2 * n - 1 - M) * h)

    )

def B(n, m, f, T, M, h):

    return (

        np.cos(2 * np.pi * f * T - alpha(n, m, M, h)) - phi_h(h, M) * np.cos(alpha(n, m, M, h))

    ) / (

        1

        + pow(phi_h(h, M), 2)

        - 2 * phi_h(h, M) * np.cos(2 * np.pi * f * T)

    )

def s_v(f, T, M, h):

    sum_a = 0

    sum_ab = 0

    for n in range(1, M + 1):