四根发射天线、一根接收天线,QPSK,4状态空时格码,一条BER曲线从0 dB扫到12 dB,你大概要跑多久?我第一次把这个4天线STTC仿真MATLAB工程搭起来的时候,答案是一整夜,而且第二天早上看到的曲线还是一条几乎不下降的斜线。问题不在信道,也不在信噪比定义,而是维特比译码里的状态回溯写反了。后来我把编码器、信道、译码器和统计脚本拆成四个独立模块,先用无噪声回环验证,再逐段加噪,BER/FER曲线才终于跟理论斜率对上。这篇就按我当时调试的顺序,把4天线STTC从码表、QPSK映射、准静态瑞利信道、维特比译码到BER/FER统计的完整工程拆开讲。适合已经写过基础MIMO仿真、但一加译码器就卡住的人,也适合想直接拿一套可复现MATLAB框架去改码表、改天线数的读者。里面会给出我实际用的参数、代码骨架、复杂度估算和几条踩坑记录,照着搭能少走不少弯路。
1. 先把链路搭对:4天线STTC仿真到底在算什么
1.1 为什么我坚持用STTC而不是拿STBC凑数
空时格码(Space-Time Trellis Code,STTC)和正交空时分组码(OSTBC)经常被放在一起比较。如果你的目标只是验证四天线分集阶数,OSTBC确实更省事,编码就是矩阵共轭转置,译码只做线性合并,十几行代码就能跑出曲线。但四发射天线的复正交设计存在码率损失,常见的4发1收OSTBC码率只有3/4,也就是说同样带宽下,每个信道用只传3个符号而不是4个。STTC是速率1的格码,每个符号周期四个天线各发一个QPSK符号,频谱效率上没有这块短板。
更关键的是,STTC把编码和分集绑在一起做联合设计。它不像OSTBC那样只追求正交性,而是在格图上选择一条条码字路径,让任意两条合法路径之间的差异矩阵尽量满秩,从而同时拿到分集增益和编码增益。对于四发一收,理论分集阶数是4,BER曲线在高信噪比段的斜率应该接近4。你如果只做OSTBC,也能看到斜率4,但曲线位置会偏右,因为缺少编码增益。STTC的维特比译码复杂度更高,但正好适合拿来练手译码器、度量计算和路径回溯,这也是我选择它做工程解析的原因。
还有一个现实考虑:很多论文和教材在讲空时码时,STTC的译码部分是最容易被一笔带过的。编码端给个生成矩阵,接收端说一句“用维特比译码”就结束了。真正动手的人会卡在状态索引、分支度量、回溯起点这些细节上。我把这条链路完整跑通一遍,后面换成8状态、16状态码,或者换成两根接收天线,框架基本不用动,只需要替换码表和改维度。这种可扩展性比单纯跑一条OSTBC曲线更有价值。
1.2 一条最小可用链路的模块清单
我把整个仿真拆成七块:比特源、符号映射、STTC编码、功率归一化、准静态瑞利信道、AWGN、维特比译码,最后是BER/FER统计。每一块的输入输出维度必须提前约定清楚,否则后面调试时会出现“信道矩阵转置了但曲线看起来只是差一点”的隐蔽错误。我的约定是:一帧包含128个信息符号加1个尾符号,每个符号携带2比特,所以每帧256个信息比特;发射天线数Nt=4,接收天线数Nr=1;编码器输出是4×130的复矩阵,每一列是一个时刻四个天线的发送符号。
| 模块 | 输入维度 | 输出维度 | 关键约定 |
|---|---|---|---|
| 比特源 | 1×256 | 1×256 | 0/1随机,按帧独立生成 |
| 符号映射 | 1×256 | 1×128 | 每2比特映射为一个QPSK符号索引0~3 |
| STTC编码 | 1×129 | 4×129 | 含1个尾符号,状态从0开始回到0 |
| 功率归一化 | 4×129 | 4×129 | 每个天线除以sqrt(Nt),总功率归一 |
| 瑞利信道 | 4×129 | 1×129 | 准静态,一帧内H不变 |
| AWGN | 1×129 | 1×129 | 复噪声,实部虚部各N0/2 |
| 维特比译码 | 1×129 + H | 1×129 | 输出符号索引,再拆成比特 |
| 统计 | 比特/帧 | BER/FER | 按帧累计错误 |
这张表看起来简单,但每一项都对应一个容易出错的点。比如“每2比特映射为一个符号索引”这一步,如果你用格雷映射,符号索引和比特的对应关系就不是简单的二进制拼接,后面算BER时又得反向查表。我在这个工程里故意用了自然映射,让符号索引的模4运算直接对应QPSK的相位旋转,编码器的线性性质更好验证。格雷映射的BER会略好一点,但码表系数需要重新校验,不能直接套。
1.3 工程目录与数据流约定
我习惯把每个功能拆成独立函数,主脚本只负责循环和统计。目录结构大概是这样:
sttc4x1/ main_ber.m sttc_encode.m sttc_decode.m build_trellis.m qpsk_map.m gen_channel.m utils/ bits2sym.m sym2bits.m主脚本main_ber.m负责SNR扫描、帧循环、误码累计和画图;sttc_encode.m只做编码,输入符号索引序列,输出四天线复符号矩阵;sttc_decode.m只做维特比译码,输入接收信号和信道矩阵,输出符号索引估计;build_trellis.m把码表展开成16条分支的输出向量和下一状态,供译码器查表。这样拆的好处是,编码器和译码器可以单独做无噪声回环测试,不需要每次跑完整的BER循环。
数据流上我坚持用“符号索引”作为中间表示,而不是直接传复数。原因是编码器的模4运算、状态转移、译码器的分支索引都建立在0~3的整数上,复数只在映射之后出现。如果你在编码器内部就把索引转成复数,后面做模4加法还得先转回来,容易把相位和索引搅混。我的做法是:bits2sym把比特变成0~3的符号索引,sttc_encode输出0~3的编码符号索引,qpsk_map再统一转成复符号并做功率归一化。译码器输出的也是0~3的索引,最后用sym2bits还原比特。整条链路只有两个地方出现复数:信道和映射之后。
2. 编码端:4状态STTC的码表、映射与功率归一化
2.1 生成矩阵与移位寄存器的对应关系
四状态STTC的状态数等于4,意味着编码器记忆的是一个QPSK符号,也就是2比特。当前输入符号记为u(t),上一时刻输入符号记为u(t-1),状态就是u(t-1)的四种取值。每个天线上的输出符号由当前输入和上一状态共同决定,写成模4运算就是:
x_i(t) = ( u(t) + g_i * u(t-1) ) mod 4,i = 1,2,3,4
其中g_i是第i个天线对应的生成系数,取0~3的整数。这个式子本质上是把输入符号和延迟一个时刻的输入符号做相位旋转再加权,四个天线用不同的g_i产生不同的码字分量。我的测试工程里用了一组简单的系数:
| 天线 | g_i | 输出表达式 |
|---|---|---|
| 1 | 1 | x1 = u(t) + u(t-1) mod 4 |
| 2 | 2 | x2 = u(t) + 2u(t-1) mod 4 |
| 3 | 3 | x3 = u(t) + 3u(t-1) mod 4 |
| 4 | 0 | x4 = u(t) |
这组系数不是为了追求最优编码增益,而是为了让四条支路的输出差异足够明显,方便我在调试时观察状态转移。换成经典文献里的四状态码表时,只需要把g_i替换掉,编码器和译码器的框架完全不用改。这里要提醒一句:如果你直接把从别处抄来的生成矩阵塞进这个式子,先确认对方的符号映射是自然映射还是格雷映射,模4运算的定义是否一致。映射不一致,码字的距离特性就变了,曲线可能比预期差很多。
状态转移也很简单:下一状态等于当前输入符号,因为编码器只记忆一个符号。四状态对应0、1、2、3四个值,每个状态有4条输入分支,总共16条分支。这个数字很小,译码器循环写起来不费劲,但也正因为小,状态索引写错一位就很容易被忽略,曲线只是“看起来差一点”,不会直接报错。
2.2 QPSK符号映射与平均功率校准
QPSK映射我采用自然相位顺序:符号索引0对应(1+1j)/sqrt(2),1对应(-1+1j)/sqrt(2),2对应(-1-1j)/sqrt(2),3对应(1-1j)/sqrt(2)。四个点都在单位圆上,相邻索引相差90度。这样符号索引的模4加法就等价于相位旋转,编码器的线性性质可以直接用相位关系验证。比特到符号的映射也用自然二进制:00→0,01→1,10→2,11→3,也就是u = 2*b1 + b2。这样一来,编码、映射、比特还原三处的索引定义完全一致,调试时少一层心智负担。
功率归一化是另一个容易埋雷的地方。四根天线同时发射,如果每个天线都发单位功率符号,总发射功率就是4。为了让不同天线数之间的BER曲线可以横向比较,我统一按总功率归一化:每个天线符号再除以sqrt(Nt)=2,保证E[||x||^2]=1。此时噪声方差N0直接由SNR定义给出,SNR_dB = 10*log10(1/N0),接收端每根天线的总信噪比就是SNR。这个约定下,四发一收和单发一收在相同SNR下的总发射功率一致,曲线才有可比性。
如果你习惯让每个天线发单位功率,那也完全可以,但要在噪声方差或SNR定义上同步调整,否则曲线会整体偏移6 dB。我见过不少人把“每天线单位功率”和“总功率归一化”混着用,最后发现译码增益比预期少6 dB,查了一整天以为是译码器问题。实际就是功率口径不统一。我的建议是:在脚本最上面写一行注释,明确“总功率归一化,每天线除以sqrt(Nt)”,后面所有噪声和SNR计算都围绕这一行来。
2.3 帧结构、状态复位与尾比特处理
帧长我设为130个符号,其中前128个是信息符号,最后1个是尾符号。尾符号固定为0,作用是让编码器状态从0出发,经过128个信息符号后,再经过一个0输入回到状态0。四状态编码器只记忆一个符号,所以一个尾符号就够了。译码时从已知的终止状态0开始回溯,不需要在四个状态里挑路径度量最小的那个,逻辑更干净。
尾符号的开销是1/130,约0.77%,对BER曲线影响很小,但对FER统计有实际意义。FER的定义是“一帧中只要有任意一个比特错误,就算一帧错误”。如果没有尾符号,译码器必须在四个终止状态里选一个度量最小的,选错终止状态会导致整帧后半段回溯全错,FER会偏高。加上尾符号之后,终止状态确定,回溯起点唯一,FER统计更稳定。对于更长的帧,比如1024个符号,尾符号开销可以忽略,但四状态码的尾符号仍然只需要1个,因为状态数只有4。
这里还有一个细节:尾符号本身也参与编码和传输,译码器会把它译出来,但统计BER时不计入尾符号对应的比特。我的做法是只比较前128个信息符号对应的256个比特,尾符号译错不影响BER,但会影响FER吗?严格来说,如果尾符号译错但信息比特全对,这一帧的信息是正确的,不应该算帧错误。所以FER统计也只检查信息比特。这一点在写统计循环时要明确,否则高信噪比下会出现“BER为0但FER不为0”的怪现象,其实是尾符号在捣乱。
3. 信道与接收端:维特比译码的度量怎么写才不跑偏
3.1 准静态瑞利信道的生成与CSI假设
STTC的码字跨多个符号周期,译码器需要利用多个时刻的接收信号做联合判决,所以信道模型通常选准静态瑞利衰落:一帧之内信道矩阵H保持不变,帧与帧之间独立变化。这个假设对应慢速移动场景,也让“一帧一个H”的实现变得简单。H是Nr×Nt的复高斯矩阵,每个元素独立同分布,实部虚部都是均值为0、方差为1/2的高斯变量,这样每个元素的平均功率是1。MATLAB里写成:
H = (randn(Nr,Nt) + 1j*randn(Nr,Nt)) / sqrt(2);接收信号模型是y = H*x + n,其中x是Nt×1的发送符号向量,n是Nr×1的复AWGN向量,每个元素的实部虚部方差为N0/2。译码器假设完美信道状态信息(CSI),也就是H已知。这个假设在仿真里很常见,因为STTC的译码度量需要H来计算欧氏距离;如果要做信道估计,那就是另一个话题了,会引入估计误差和导频开销,曲线会再偏移。
准静态假设下,四发一收的理论分集阶数是Nt×Nr=4。也就是说,高信噪比段BER曲线的斜率应该是4,每增加10 dB,BER大约下降四个数量级。你如果看到斜率只有2,先别怀疑码表,查一下H是不是每帧都重新生成、有没有误把每符号都变信道、或者接收天线数是不是设成了1但维度写错。斜率不对通常不是译码器算法的问题,而是信道或分集路径没对上。
3.2 分支度量的矩阵化计算
维特比译码的核心是分支度量。对于某一个时刻t,给定上一状态s和当前输入u,编码器输出一个Nt×1的符号向量x(s,u),接收端计算这个假设下的预测接收信号H*x(s,u),再和实际接收信号y比较,取欧氏距离的平方:
M(s,u) = || y(t) - H * x(s,u) ||^2
四状态、四输入,一共16条分支。每条分支算一次Nr×Nt的矩阵向量乘,如果放在三重循环里,MATLAB的解释器开销会很大。我的做法是预先把16条分支的输出向量拼成一个Nt×16的矩阵X,然后一次矩阵乘法算出所有预测:
X = zeros(Nt, 16); for s = 0:3 for u = 0:3 idx = s*4 + u + 1; x_idx = mod(u + G*s, 4); X(:,idx) = qpsk_map(x_idx) / sqrt(Nt); end end % 每个时刻一次矩阵乘法 pred = H * X; % Nr x 16 metric = sum(abs(y_t - pred).^2, 1); % 1 x 16这段代码是整个译码器提速的关键。H是Nr×Nt,X是Nt×16,一次乘法得到Nr×16的预测矩阵,再用隐式扩展减去y_t,取模平方,沿第一维求和,得到1×16的度量向量。16条分支的度量一次算完,避免了在状态循环里反复调用矩阵乘法。对于Nr=1、Nt=4的情况,单次乘法只有64次复数乘加,130个符号一帧也就八千多次,MATLAB跑起来很快。
这里要注意metric的排列顺序。我按idx = s*4 + u + 1排列,s是上一状态,u是当前输入。译码器在ACS循环里取度量时,必须用同样的索引方式,否则会出现“度量算对了但分支取错”的问题。这种错误不会报错,只会让曲线变差,非常隐蔽。我的习惯是在build_trellis.m里同时返回一个16×2的branch表,第一列是上一状态,第二列是当前输入,译码器直接按表取,不手算索引。
3.3 加比选与回溯:状态索引和幸存路径的坑
加比选(ACS)的逻辑不复杂:对于每个下一状态,遍历所有可能到达它的上一状态,把上一状态的路径度量加上分支度量,取最小的那条作为幸存路径,同时记录幸存的前一状态和当前输入。四状态码里,下一状态等于当前输入u,所以每个下一状态有4条入分支。我的实现是外层循环上一状态s,内层循环输入u,算出下一状态ns=u,然后更新ns的路径度量。这样写的好处是分支表和编码器的状态转移完全一致,不容易绕晕。
pm = [0; inf; inf; inf]; % 初始状态为0 surv_s = zeros(L, 4); % 每个时刻、每个下一状态的幸存前一状态 surv_u = zeros(L, 4); % 每个时刻、每个下一状态的幸存输入 for t = 1:L pm_new = inf(4,1); for s = 0:3 for u = 0:3 ns = u; % 下一状态等于当前输入 idx = s*4 + u + 1; m = metric(t, idx); if pm(s+1) + m < pm_new(ns+1) pm_new(ns+1) = pm(s+1) + m; surv_s(t, ns+1) = s; surv_u(t, ns+1) = u; end end end pm = pm_new; end回溯从终止状态0开始,倒着读surv_u:
state = 0; u_hat = zeros(1, L); for t = L:-1:1 u_hat(t) = surv_u(t, state+1); state = surv_s(t, state+1); end这里有两个坑。第一个是MATLAB的1-based索引:状态0对应数组下标1,状态3对应下标4。如果你在存surv_s时用了s而不是s+1,或者回溯时用了state而不是state+1,就会取到错位的幸存路径。第二个是回溯方向:surv_u(t, ns+1)记录的是“在时刻t到达下一状态ns的那条分支的输入”,回溯时从最后一刻的终止状态开始,先用surv_u取出输入,再用surv_s跳到前一状态。顺序反了,整条路径就乱了。我第一次跑出那条不下降的曲线,就是这两个坑叠在一起。
3.4 复杂度估算与截断深度选择
四状态STTC的译码复杂度很低,但把数字算清楚有助于你判断跑一条曲线需要多久。每个符号时刻16条分支,每条分支计算Nr×Nt的欧氏距离,Nr=1、Nt=4时是4次复数乘加,16条分支就是64次。一帧130个符号约8320次,一万帧就是8320万次复数运算。MATLAB的矩阵化实现下,这个量级在普通笔记本上跑一条12个SNR点的曲线大概几十分钟到一小时,具体取决于帧数和错误目标。
| 项目 | 数值 | 说明 |
|---|---|---|
| 状态数 | 4 | 记忆一个QPSK符号 |
| 每状态输入分支 | 4 | QPSK有4个符号 |
| 每符号总分支 | 16 | 4状态×4输入 |
| 每分支复数乘加 | Nr×Nt | Nr=1,Nt=4时为4 |
| 每符号复数乘加 | 64 | 16×4 |
| 每帧(130符号) | 8320 | 不含ACS比较 |
| 一万帧 | 8320万 | 粗略量级 |
回溯深度方面,四状态码的约束长度是2,理论上回溯深度取约束长度的5倍左右就够了,也就是10个符号。我用15个符号作为截断深度,在长帧仿真里完全够用。但这个工程用的是130符号的短帧,直接全帧回溯更简单,也不需要担心截断带来的性能损失。如果你把帧长加到几千个符号,建议改成滑动窗口回溯,每15个符号输出一次判决,内存占用和计算量都会更友好。
4. MATLAB实操:从主循环到BER/FER统计
4.1 参数区与SNR扫描脚本
主脚本我习惯把所有可调参数放在最上面,避免在循环里出现魔法数字。下面是我实际用的骨架,SNR从0到12 dB,步长2 dB,每个SNR点最多跑20000帧,或者累计到100个帧错误就停。低信噪比下错误多,很快就停;高信噪比下错误少,跑满20000帧,曲线的高端会平一些,但至少不会出现“零错误”的假点。
clear; clc; rng(2024); Nt = 4; Nr = 1; M = 4; bitsPerSym = 2; L_info = 128; L = L_info + 1; % 加1个尾符号 SNR_dB = 0:2:12; maxFrames = 20000; targetFrameErr = 100; G = [1; 2; 3; 0]; % 四天线生成系数 qpsk = [1+1j, -1+1j, -1-1j, 1-1j] / sqrt(2); trellis = build_trellis(G, qpsk, Nt); ber = zeros(size(SNR_dB)); fer = zeros(size(SNR_dB)); for iSNR = 1:length(SNR_dB) snr = SNR_dB(iSNR); N0 = 10^(-snr/10); frameErr = 0; bitErr = 0; totalBits = 0; totalFrames = 0; while totalFrames < maxFrames && frameErr < targetFrameErr bits = randi([0 1], 1, L_info*bitsPerSym); u = bits2sym(bits); % 1 x L_info u = [u, 0]; % 尾符号 x = sttc_encode(u, G, qpsk, Nt); % Nt x L H = (randn(Nr,Nt) + 1j*randn(Nr,Nt)) / sqrt(2); n = sqrt(N0/2) * (randn(Nr,L) + 1j*randn(Nr,L)); y = H * x + n; u_hat = sttc_decode(y, H, trellis, L, Nt); bits_hat = sym2bits(u_hat(1:L_info)); err = sum(bits ~= bits_hat); bitErr = bitErr + err; if err > 0 frameErr = frameErr + 1; end totalBits = totalBits + length(bits); totalFrames = totalFrames + 1; end ber(iSNR) = bitErr / totalBits; fer(iSNR) = frameErr / totalFrames; fprintf('SNR=%2d dB, frames=%6d, BER=%.3e, FER=%.3e\n', ... snr, totalFrames, ber(iSNR), fer(iSNR)); end semilogy(SNR_dB, ber, 'o-', SNR_dB, fer, 's--'); grid on; xlabel('SNR (dB)'); ylabel('BER / FER'); legend('BER', 'FER');这段脚本里bits2sym和sym2bits是互逆的查表函数。自然映射下,u = 2*b1 + b2,反变换就是b1 = floor(u/2),b2 = mod(u,2)。不要小看这两个函数,如果它们的映射方向和QPSK映射表不一致,BER会莫名其妙地变差,但FER可能看起来正常,因为帧错误主要由符号错误引起,而比特错误的分布被映射关系扭曲了。
4.2 编码函数实现
编码函数只做一件事:按状态转移生成四天线符号索引,再映射成复符号。输入是1×L的符号索引,输出是Nt×L的复矩阵。状态从0开始,每个时刻用当前状态(上一输入)和当前输入算四个输出,然后状态更新为当前输入。
function x = sttc_encode(u, G, qpsk, Nt) L = length(u); x = zeros(Nt, L); s = 0; % 初始状态 for t = 1:L for i = 1:Nt x_idx = mod(u(t) + G(i)*s, 4); x(i,t) = qpsk(x_idx + 1) / sqrt(Nt); end s = u(t); end end这个函数很短,但有一个细节:功率归一化放在映射之后、编码输出之前,除以sqrt(Nt)。这样编码器输出的每一列总功率为1,后面的信道和噪声计算可以直接用N0。如果你把归一化放到信道之后,或者忘了除以sqrt(Nt),SNR的口径就变了。我建议把归一化写死在编码函数里,主脚本不再重复处理,避免两处都归一化导致功率被压两次。
4.3 译码函数实现
译码函数是工程的核心。输入是Nr×L的接收信号、Nr×Nt的信道矩阵、预构建的trellis结构,输出1×L的符号索引估计。trellis里存了16条分支的输出向量、下一状态和分支索引,译码循环直接查表。
function u_hat = sttc_decode(y, H, trellis, L, Nt) Nr = size(H,1); ns = 4; pm = [0; inf; inf; inf]; surv_s = zeros(L, ns); surv_u = zeros(L, ns); metric = zeros(1, 16); % 预构建的16条分支输出 X = trellis.X; % Nt x 16 for t = 1:L pred = H * X; % Nr x 16 metric = sum(abs(y(:,t) - pred).^2, 1); pm_new = inf(ns,1); for s = 0:3 for u = 0:3 idx = s*4 + u + 1; ns_next = trellis.next(idx); m = metric(idx); if pm(s+1) + m < pm_new(ns_next+1) pm_new(ns_next+1) = pm(s+1) + m; surv_s(t, ns_next+1) = s; surv_u(t, ns_next+1) = u