一、正文

在编译时出现了一个小问题，我的 .M 文件无法正常打开，显示是不受支持的 MATLAB 代码，后来我重新定义为小写的 .m 文件，问题得到了解决。

当然，另一种解决方案是将函数内容直接写道主函数里面，此处不表。

k = menu('pi/4 QPSK Plot Options',...
        'Unfiltered pi/4 QPSK Signal Constellation',...
        'Unfiltered pi/4 QPSK Eye Diagram',...
        'Filtered pi/4 QPSK Signal Constellation',...
        'Filtered pi/4 OQPSK Eye Diagram',...
        'Unfiltered Direct and Quadrature Signals',...
        'Filtered Direct and Quadrature Signals',...
        'Exit Program');

这段代码创建了一个菜单，其中包含了7个选项，每个选项对应一个图形的类型。用户可以通过选择不同的选项来查看不同的图形。如果选择“Exit Program”，则程序将退出。

menu函数参数表中的 “…” 有什么用？

在这里，“…” 的作用是让你可以在同一个 menu 函数调用中提供多个选项，而不需要在同一行中输入所有选项。

%
% set direct and quadrature bit streams
%
data = round(rand(1,bits));
dd = data(1:2:bits-1);
qq = data(2:2:bits);

串井转换器交替将（奇序数的）符号分配给同相信道，将其他（偶序数的）符号分配给正交信道，千是有

a
(
1
)
a
(
2
)
a
(
3
)
a
(
4
)
⋯
a
(
k
)
⋯
=
d
(
1
)
q
(
1
)
d
(
2
)
q
(
2
)
⋯
d
(
k
+
1
2
)
q
(
k
2
+
1
)
⋯
a(1) a(2) a(3) a(4) \cdots a(k) \cdots=d(1) q(1) d(2) q(2) \cdots d\left(\frac{k+1}{2}\right) q\left(\frac{k}{2}+1\right) \cdots
a(1)a(2)a(3)a(4)⋯a(k)⋯=d(1)q(1)d(2)q(2)⋯d(2k+1​)q(2k​+1)⋯

这段代码是用随机数生成一个长度为bits的二进制数据流，并将其分成了直接通道（dd）和象限通道（qq）两个部分。其中dd包含了数据流中的奇数位，qq包含了数据流中的偶数位。这个过程是为了将二进制数据流转换成复数信号，以便进行QPSK调制。

theta(1) = iphase;						% set initial phase
thetaout(1:sps) = theta(1)*ones(1,sps);

$$\theta (k)=\theta (k-1)+\varphi (k)$$

for k=2:m
    if dd(k) == 1
        phi_k = (2*qq(k)-1)*pi/4;
    else
        phi_k = (2*qq(k)-1)*3*pi/4;
    end
    theta(k) = phi_k + theta(k-1);
    for i=1:sps
        j = (k-1)*sps+i;
        thetaout(j) = theta(k);
    end
end
d = cos(thetaout);
q = sin(thetaout);

这一段代码是有问题的，因为教材P195 说的很清楚，且初始相位定义为 $\theta (0)=0$。按照上面 MATLAB 代码的写法，dd(1)和qq(1)是永远用不到的。我们知道就好，不影响整体程序的运行“应该”。

1. 星座图

% File: sigcon.m
% Software given here is to accompany the textbook: W.H. Tranter,
% K.S. Shanmugan, T.S. Rappaport, and K.S. Kosbar, Principles of
% Communication Systems Simulation with Wireless Applications,
% Prentice Hall PTR, 2004.
% 
function []=sigcon(x,y)
plot(x,y);
axis('square')
axis('equal')
xlabel('Direct Channel')
ylabel('Quadrature Channel')
end
% End of function file.

x表示同相信号，y表示正交信号。

2. 眼图

% File: dqeye.m
% Software given here is to accompany the textbook: W.H. Tranter,
% K.S. Shanmugan, T.S. Rappaport, and K.S. Kosbar, Principles of
% Communication Systems Simulation with Wireless Applications,
% Prentice Hall PTR, 2004.
%
% plot unfiltered eye diagram
%xd是同相信号，xq是正交信号，m在主程序中是4*sps=40
function [] = dqeye(xd,xq,m)
lx = length(xd);                % samples in data segment，选取200个symbol,每个symbol取10个点，lx=2000
kcol = floor(lx/m);             % number of columns
xda = [0,xd]; xqa = [0,xq];     % append zeros,默认初始相位是0
for j = 1:kcol                  % column index
    for i = 1:(m+1)             % row index
        kk = (j-1)*m+i;         % sample index
        y1(i,j) = xda(kk);
        y2(i,j) = xqa(kk);
    end
end
subplot(211)                    % direct channel
plot(y1);
title('D/Q EYE DIAGRAM');
xlabel('Sample Index');
ylabel('Direct');
subplot(212)                    % quadrature channel
plot(y2);
xlabel('Sample Index');
ylabel('Quadratute');
subplot(111)                    % restore
% End of function file.

对plot(y2); 的一点解释:

在Matlab中，当使用plot(a)绘制一个矩阵时，Matlab会自动将该矩阵的每一列视为一个向量，并将其绘制在一个单独的坐标系上，横坐标默认为该列在矩阵中的 row 索引，纵坐标则为该列的值。因此，对于一个40*50的矩阵，Matlab会绘制出50条线段，每条线段有40个点，横坐标分别为1到40。

眼图这里的代码有点难以理解，我来做一点解释。

d = cos(thetaout);
q = sin(thetaout);
elseif k ==2
        dqeye(d,q,4*sps)        	% plot unfiltered eye diagram
        pause

我们通过上面这个代码得到的d和q实际上是 $1\times 2000$ 的向量，选取200个symbol,每个symbol取10个点，lx=2000。那上图中一个 Segment 在程序中我们取 4*sps，如下图示意。

但是如果真的取4*sps画图的话，因为我们这里是模拟仿真，会发现眼图是恰好会重合在一起，为了示意眼图，我们用下面的方式排列矩阵并且绘图，就能得到正确的眼图了。

8.2.3 增益、延迟和信噪比

实低通信号的理论推导

对线性时不变无失真系统，系统中任意点信号 $y(t)$ 是对输入参考信号 $x(t)$ 按幅度比例缩放并加入延迟后所得的信号。因此， 我们可以将无失真信号写为：

$$z(t)=Ax(t-\tau )$$

其中， $A$ 是增益， $\tau$ 是系统中定义SNR的那个点的群时延。令 $x(t)$ 为参考信号， $y(t)$ 为观测信号，于是
$$y(t)=Ax(t-\tau )+n(t)+d(t)\text{\hspace{1em}(21)}$$

其中， $n(t)$ 代表外部加性噪声， $d(t)$ 是系统引人的跟信号相关的内部失真， 这种内部失真可能由码间千扰或非线性造成。描述系统中不同信号之间关系的方框图如图8-11 所示。

噪声功率定义为 $y(t)$ 和无失真系统输出 $z(t)=Ax(t-\tau )$之间的MSE , 即：

ε
(
A
,
τ
)
=
E
{
[
y
(
t
)
−
A
x
(
t
−
τ
)
]
2
}
\varepsilon(A, \tau)=E\left\{[y(t)-A x(t-\tau)]^{2}\right\}
ε(A,τ)=E{[y(t)−Ax(t−τ)]2}

$A$ 和 $\tau$ 理想估计值是最小化 $\epsilon (A,\tau )$的值。上述表达式可写为：

ε
(
A
,
τ
)
=
E
{
y
2
(
t
)
+
A
2
x
2
(
t
−
τ
)
−
2
A
x
(
t
−
τ
)
y
(
t
)
}
(23)
\varepsilon(A, \tau)=E\left\{y^{2}(t)+A^{2} x^{2}(t-\tau)-2 A x(t-\tau) y(t)\right\}\tag{23}
ε(A,τ)=E{y2(t)+A2x2(t−τ)−2Ax(t−τ)y(t)}(23)

对于平稳信号， 各阶矩独立于时间起点。此外， 和的期望等于期望的和，因此， 方程(8-23)可写为：

ε
(
A
,
τ
)
=
E
{
y
2
(
t
)
}
+
A
2
E
{
x
2
(
t
)
}
−
2
A
E
{
x
(
t
)
y
(
t
+
τ
)
}
(24)
\varepsilon(A, \tau)=E\left\{y^{2}(t)\right\}+A^{2} E\left\{x^{2}(t)\right\}-2 A E\{x(t) y(t+\tau)\}\tag{24}
ε(A,τ)=E{y2(t)}+A2E{x2(t)}−2AE{x(t)y(t+τ)}(24)

这里的 $E\{x(t)y(t+\tau )\}$ 是求得最大的期望相关功率，在仿真时我们选择近似用 $x(t)$ 的长度来代替时长。

因此，信噪比为

S
N
=
R
x
y
2
(
τ
m
)
P
x
[
P
x
P
x
P
y
−
R
x
y
2
(
τ
m
)
]
=
R
x
y
2
(
τ
m
)
P
x
P
y
−
R
x
y
2
(
τ
m
)
(30)
\frac{S}{N}=\frac{R_{x y}^{2}\left(\tau_{m}\right)}{P_{x}}\left[\frac{P_{x}}{P_{x} P_{y}-R_{x y}^{2}\left(\tau_{m}\right)}\right]=\frac{R_{x y}^{2}\left(\tau_{m}\right)}{P_{x} P_{y}-R_{x y}^{2}\left(\tau_{m}\right)}\tag{30}
NS​=Px​Rxy2​(τm​)​[Px​Py​−Rxy2​(τm​)Px​​]=Px​Py​−Rxy2​(τm​)Rxy2​(τm​)​(30)

$x(t)$ 和 $y(t)$ 的相关系数 $\rho$ 定义为

ρ
=
R
x
y
(
τ
m
)
P
x
P
y
(31)
\rho=\frac{R_{x y}\left(\tau_{m}\right)}{\sqrt{P_{x} P_{y}}}\tag{31}
ρ=Px​Py​​Rxy​(τm​)​(31)

根据这个定义，在测量点 $y(t)$ 处， SNR的表达式很简单：

S
N
=
ρ
2
1
−
ρ
2
(32)
\frac{S}{N}=\frac{\rho^{2}}{1-\rho^{2}}\tag{32}
NS​=1−ρ2ρ2​(32)

我们都知道圆周卷积有这个公式

y
(
n
)
=
[
∑
m
=
L
−
1
x
1
(
m
)
∑
r
=
−
∞
∞
x
2
(
n
+
r
L
−
m
)
]
R
L
(
n
)
=
[
∑
r
=
−
∞
∞
∑
m
=
L
−
1
x
1
(
m
)
x
2
(
n
+
r
L
−
m
)
]
R
L
(
n
)
=
[
∑
r
=
−
∞
∞
y
l
(
n
+
r
L
)
]
R
L
(
n
)
\begin{aligned} y(n) & =\left[\sum_{m=0}^{L-1} x_{1}(m) \sum_{r=-\infty}^{\infty} x_{2}(n+r L-m)\right] \mathrm{R}_{\mathrm{L}}(n) \\ & =\left[\sum_{r=-\infty}^{\infty} \sum_{m=0}^{L-1} x_{1}(m) x_{2}(n+r L-m)\right] \mathrm{R}_{L}(n)\\ &=\left[\sum_{r=-\infty}^{\infty} y_{l}(n+r L)\right] R_{L}(n) \end{aligned}
y(n)​=[m=0∑L−1​x1​(m)r=−∞∑∞​x2​(n+rL−m)]RL​(n)=[r=−∞∑∞​m=0∑L−1​x1​(m)x2​(n+rL−m)]RL​(n)=[r=−∞∑∞​yl​(n+rL)]RL​(n)​
我们知道圆周卷积有定义

r
ˉ
x
y
(
m
)
=
∑
n
=
N
−
1
x
(
n
)
y
(
(
n
−
m
)
)
N
R
N
(
m
)
=
∑
n
=
N
−
1
x
(
n
)
y
(
(
−
(
m
−
n
)
)
)
N
R
N
(
m
)
=
x
(
m
)
N
◯
y
(
N
−
m
)
\begin{aligned} \bar{r}_{x y}(m) & =\sum_{n=0}^{N-1} x(n) y((n-m))_{N} R_{N}(m) \\ & =\sum_{n=0}^{N-1} x(n) y((-(m-n)))_{N} R_{N}(m)=x(m) \textcircled{N} y(N-m) \end{aligned}
rˉxy​(m)​=n=0∑N−1​x(n)y((n−m))N​RN​(m)=n=0∑N−1​x(n)y((−(m−n)))N​RN​(m)=x(m)N◯y(N−m)​
那么必然圆周卷积也是线性相关的平移，即

r
ˉ
x
y
(
m
)
=
∑
n
=
N
1
−
1
x
(
n
)
y
(
(
n
−
m
)
)
N
R
N
(
m
)
=
[
∑
n
=
N
1
−
1
x
(
n
)
∑
r
=
−
∞
∞
y
(
n
−
m
+
r
L
)
]
R
N
(
m
)
=
[
∑
r
=
−
∞
∞
∑
n
=
n
=
N
1
−
1
x
(
n
)
y
(
n
−
(
m
−
r
L
)
)
]
R
N
(
m
)
=
[
∑
r
=
−
∞
∞
r
x
y
(
n
−
r
N
)
]
R
N
(
m
)
\begin{aligned} \bar{r}_{x y}(m) & =\sum_{n=0}^{N_1-1} x(n) y((n-m))_{N} R_{N}(m) \\ &=[\sum_{n=0}^{N_1-1} x(n) \sum_{r=-\infty}^{\infty}y(n-m+rL) ]R_{N}(m)\\ &=[\sum_{r=-\infty}^{\infty} \sum_{n=0}^{n=N_1-1}x(n)y(n-(m-rL))]R_{N}(m)\\ &=[\sum_{r=-\infty}^{\infty}r_{xy}(n-r N)] R_{N}(m) \end{aligned}
rˉxy​(m)​=n=0∑N1​−1​x(n)y((n−m))N​RN​(m)=[n=0∑N1​−1​x(n)r=−∞∑∞​y(n−m+rL)]RN​(m)=[r=−∞∑∞​n=0∑n=N1​−1​x(n)y(n−(m−rL))]RN​(m)=[r=−∞∑∞​rxy​(n−rN)]RN​(m)​
其中，
r
x
y
(
n
)
=
x
(
n
)
∗
y
(
−
n
)
r_{xy}(n)=x(n)*y(-n)
rxy​(n)=x(n)∗y(−n) 是线性相关。

% File: snrmse.m
% Software given here is to accompany the textbook: W.H. Tranter, 
% K.S. Shanmugan, T.S. Rappaport, and K.S. Kosbar, Principles of 
% Communication Systems Simulation with Wireless Applications, 
% Prentice Hall PTR, 2004.
%
function [gain,delay,px,py,rxy,rho,snrdb] = snrmse(x,y)
ln = length(x);                 % Set length of the reference (x) vector
fx = fft(x,ln);                 % FFT the reference (x) vector
fy = fft(y,ln);                 % FFT the measurement (y) vector
fxconj = conj(fx);              % Conjugate the FFT of the reference vector
sxy = fy .* fxconj;             % Determine the cross PSD 
rxy = ifft(sxy,ln);             % Determine the cross correlation function
rxy = real(rxy)/ln;             % Take the real part and scale
px = x*x'/ln;                   % Determine power in reference vector
py = y*y'/ln;                   % Determine power in measurement vector
[rxymax,j] = max(rxy);          % Find the max of the crosscorrelation
gain = rxymax/px;               % Here's the gain
delay = j-1;                    % Here's the delay
rxy2 = rxymax*rxymax;           % Square rxymax for later use
rho = rxymax/sqrt(px*py);       % Here's the correlation coefficient
snr = rxy2/(px*py-rxy2);        % Here's the snr
snrdb = 10*log10(snr);          % Here's the snr in db
% End of script file.

这个MATLAB代码是很有问题的，我来解释一下。x 是输入参考信号 $x(t)$ ，$y(t)$ 为观测信号。

正常来说我们想要用 FFT 做线性相关必须要先补零，但是这里没有补零，那么我们直接用 ln = length(x); 的信号会得到什么？

没错，我们得到的实际上是循环卷积。

那么为什么结果依然还是基本上没有问题的呢？如下图

其实问题出在 $z(t)=Ax(t-\tau )$ 上，这里 $\tau$ 其实是非常小的，我们从下文可以知道 $\tau$ 是 64，其实是非常小的，相较于信号的长度来看。我们来看我们直接得到的信号是什么。

可以看到，$m=-64$ 的时候会取到相关系数的最大值，此时 $m$ 较大值处的影响是比较小的，可以近似的认为这就是
r
x
y
r_{xy}
rxy​ 的值。

此外还有一点，这里取实部是没有必要的，因为实函数的相关函数必然也是个实数。

%在输入x和y都是实数的情况下,相关函数rxy必定是实数，这里是没有必要的
rxy = ifft(sxy,ln);             % Determine the cross correlation function
rxy = real(rxy)/ln;             % Take the real part and scale
% isreal(rxy)
% ans =
% 
%   logical
% 
%    1
%

二、附录

r
ˉ
x
y
(
m
)
=
∑
n
=
N
1
−
1
x
(
n
)
y
(
(
n
−
m
)
)
N
R
N
(
m
)
=
[
∑
n
=
N
1
−
1
x
(
n
)
∑
r
=
−
∞
∞
y
(
n
−
m
+
r
L
)
]
R
N
(
m
)
=
[
∑
r
=
−
∞
∞
∑
n
=
n
=
N
1
−
1
x
(
n
)
y
(
n
−
(
m
−
r
L
)
)
]
R
N
(
m
)
=
[
∑
r
=
−
∞
∞
r
x
y
(
n
−
r
N
)
]
R
N
(
m
)
\begin{aligned} \bar{r}_{x y}(m) & =\sum_{n=0}^{N_1-1} x(n) y((n-m))_{N} R_{N}(m) \\ &=[\sum_{n=0}^{N_1-1} x(n) \sum_{r=-\infty}^{\infty}y(n-m+rL) ]R_{N}(m)\\ &=[\sum_{r=-\infty}^{\infty} \sum_{n=0}^{n=N_1-1}x(n)y(n-(m-rL))]R_{N}(m)\\ &=[\sum_{r=-\infty}^{\infty}r_{xy}(n-r N)] R_{N}(m) \end{aligned}
rˉxy​(m)​=n=0∑N1​−1​x(n)y((n−m))N​RN​(m)=[n=0∑N1​−1​x(n)r=−∞∑∞​y(n−m+rL)]RN​(m)=[r=−∞∑∞​n=0∑n=N1​−1​x(n)y(n−(m−rL))]RN​(m)=[r=−∞∑∞​rxy​(n−rN)]RN​(m)​

为了验算圆周相关是线性相关的平移，我花了很多时间来验证这个公式，代码分享给大家。

% 验算圆周相关
L=4;
x=[1,1,5,7];y=[2,2,3,7];
% 直接用卷积的方法求相关函数
corr=xcorr(x,y)
%用FFT的方法求相关函数
lnx = length(x);                
lny = length(y); 
fx = fft(x,lnx+lny-1);                 % FFT the reference (x) vector
fy = fft(y,lny+lnx-1);                 % FFT the measurement (y) vector
fyconj = conj(fy);             
sxy = fx .* fyconj;             
rxy = ifft(sxy,lnx+lny-1);             %计算出来的恰好以 lnx+lny-1 为长度的圆周相关
%需要将FFT得到的相关函数向右位移Ny-1个单位,得到的rxy第一个值的下标是和线性相关一致
rxy=circshift(rxy,lny-1)
rxy_1=[rxy,zeros(1,L)]%补零,为平移相加做准备
%接下来求圆周相关结果
rxy_2=[zeros(1,L),rxy(1:end)]%向右移动
rxy_3=rxy_1+rxy_2
RXY=rxy_3(lny:1:lny+L-1)%圆周相关取主值
%%
%接下来我们直接以圆周相关周期L，用FFT计算圆周相关
ln = length(L);                 % Set length of the reference (x) vector
fx = fft(x,L);                 % FFT the reference (x) vector
fy = fft(y,L);                 % FFT the measurement (y) vector
fyconj = conj(fy);              % Conjugate the FFT of the reference vector
sxy = fx .* fyconj;             % Determine the cross PSD 
rxy2 = ifft(sxy,L)             % Determine the cross correlation function