简介：这份资源面向通信工程、电子信息类专业学生及无线通信仿真入门者，提供GMSK调制解调的完整MATLAB实现方案，帮助理解GSM系统中这一高效抗干扰调制技术的原理与流程。压缩包共4个文件，均为m脚本，整体约4KB，涵盖数据生成、预处理、调制、信道模拟、解调与维特比译码等功能模块，结构紧凑便于逐段阅读与运行。资源围绕高斯滤波预滤波、π/4-QPSK符号映射、相干检测及卷积解码等关键环节展开，读者可据此搭建端到端仿真链路，观察相位平滑对多径与频偏的抑制效果，并借助维特比译码评估误码性能。目前已有2284人学习，适合作为课程设计、实验报告或自学练手的参考代码，通过修改滚降系数与信道参数，可进一步对比不同条件下的系统表现。

1. GMSK调制解调MATLAB实现：从零搭建一套能跑通的通信链路

做无线通信基带仿真的开发者，大概率绕不开GMSK这个调制方式。它频谱效率高、包络恒定、对功放非线性不敏感，GSM系统用的就是它。但真到动手写代码的时候，很多人会卡在几个地方：高斯滤波器的BT参数怎么影响眼图、Viterbi解调的状态怎么定义、相位轨迹怎么累积才不会发散。这套MATLAB实现就是冲着这些问题来的，它把GMSK的调制、信道加噪、解调完整串成一条链路，每个模块单独成函数，参数集中在脚本头部，改一个BT值就能看到误码率曲线的变化。适合正在做通信课程设计、基带算法验证，或者想从理论公式落到可运行代码的开发者。

2. 核心模块拆解：高斯滤波、相位累积与Viterbi解调

2.1 高斯低通滤波器：BT参数决定一切

GMSK的本质是在MSK调制前加一个高斯低通滤波器，用来平滑相位跳变、压缩频谱旁瓣。这个滤波器的带宽由BT参数控制，B是滤波器3dB带宽，T是码元周期。BT越小，频谱越紧凑，但码间干扰越严重，解调难度上升。常见取值是BT=0.3（GSM标准）和BT=0.5（频谱效率与误码率的折中）。

实现上，高斯滤波器的冲激响应是无限长的，实际仿真里要截断。我一般取±3个码元周期，再归一化保证直流增益为1。下面这段代码是滤波器系数生成的核心：

% 高斯低通滤波器系数生成 % BT: 带宽-码元周期积, Tb: 码元周期, L: 每码元采样点数, span: 截断符号数 function h = gaussian_filter(BT, Tb, L, span) Ts = Tb / L; % 采样周期 t = -span*Tb/2 : Ts : span*Tb/2 - Ts; % 时间向量 alpha = sqrt(2/log(2)) * pi * BT / Tb; % 高斯脉冲参数 h = alpha / sqrt(2*pi) * exp(-alpha^2 * t.^2 / 2); h = h / sum(h); % 归一化，保证增益为1 end

这里有几个参数需要说清楚。BT直接决定滤波器的陡峭程度，span控制截断长度，太小会导致频谱泄漏，太大浪费计算量。L是过采样倍数，一般取8或16，保证相位轨迹足够平滑。归一化那一步不能省，否则调制后的信号幅度会随BT变化，后续加噪的SNR就对不上了。

滤波器系数算出来之后，要对输入的NRZ（不归零）码元序列做卷积。注意输入序列要先做符号映射，把0/1映射成-1/+1，再上采样到采样率。这一步如果搞反了顺序，滤波器输出会完全不对。

2.2 相位累积与正交调制

高斯滤波后的信号还是基带脉冲，需要累积成相位才能做角度调制。GMSK的相位表达式是：

φ(t) = π/2 * ∑ a_i * ∫ g(τ - iT) dτ

其中a_i是符号，g(t)是高斯滤波后的脉冲。实际实现时，先把滤波输出做积分（累加），再乘以π/2的调制指数。这里有个容易翻车的地方：积分要用累积求和而不是MATLAB的integral函数，因为信号是离散的，用数值积分反而会引入额外误差。

% 相位累积与GMSK调制 % data: 符号序列(-1/+1), h: 高斯滤波器系数, L: 过采样倍数 function [s, phase] = gmsk_modulate(data, h, L) % 上采样并滤波 data_upsampled = upsample(data, L); filtered = conv(data_upsampled, h, 'same'); % 相位累积：先积分再乘调制指数 phase_increment = cumsum(filtered) * pi / 2; phase = phase_increment; % 正交调制 s = exp(1j * phase); end

cumsum那一步是核心，它把滤波后的脉冲变成了连续相位。pi/2是GMSK的调制指数，对应MSK的最小频移。最后取复指数得到基带信号，实部是I路，虚部是Q路。如果要做射频仿真，再乘上载波频率即可。

这里有个细节：conv的'same'选项保证输出长度和输入一致，但会引入群延迟。如果后续要做符号同步，这个延迟必须补偿。我一般会在解调端用相同的滤波器再卷一次，延迟自然抵消。

2.3 Viterbi解调：状态定义与路径度量

GMSK的解调可以用相干解调，也可以用非相干解调。这套实现用的是Viterbi最大似然序列估计，属于相干解调的一种。核心思想是把GMSK的相位状态建模成有限状态机，每个状态对应最近几个符号的组合，然后搜索使接收信号似然函数最大的路径。

状态数取决于记忆长度。BT=0.3时，高斯滤波器的影响大约持续3个符号周期，所以状态数取2^(3-1)=4或2^3=8。状态定义如下：

% Viterbi解调：状态转移与路径度量 % rx: 接收信号, h: 高斯滤波器系数, L: 过采样倍数, M: 记忆长度 function bits = gmsk_demodulate(rx, h, L, M) num_states = 2^(M-1); % 状态数 state_bits = de2bi(0:num_states-1, M-1, 'left-msb'); % 预计算每个状态的期望相位 expected_phase = zeros(num_states, L); for s = 1:num_states sym = [state_bits(s,:), 1]; % 假设下一个符号为+1 phase_tmp = cumsum(conv(upsample(sym, L), h, 'same')) * pi/2; expected_phase(s,:) = phase_tmp(end-L+1:end); end % Viterbi主循环 path_metric = zeros(num_states, 1); survivor_path = zeros(num_states, length(rx)/L); for k = 1:length(rx)/L rx_seg = rx((k-1)*L+1 : k*L); new_metric = inf(num_states, 1); new_survivor = zeros(num_states, 1); for s = 1:num_states for bit = [0, 1] % 计算分支度量 if bit == 0 exp_phase = expected_phase(s,:); else exp_phase = expected_phase(s,:) + pi/2; end branch_metric = sum(abs(rx_seg - exp(1j*exp_phase)).^2); % 状态转移 next_state = mod(s-1, 2^(M-2)) * 2 + bit + 1; total_metric = path_metric(s) + branch_metric; if total_metric < new_metric(next_state) new_metric(next_state) = total_metric; new_survivor(next_state) = s; end end end path_metric = new_metric; survivor_path(:, k) = new_survivor; end % 回溯 [~, best_state] = min(path_metric); bits = zeros(1, length(rx)/L); for k = length(rx)/L:-1:1 bits(k) = mod(best_state-1, 2); best_state = survivor_path(best_state, k); end end

这段代码是整个实现里最绕的部分，逐块解释一下。state_bits用de2bi生成所有可能的状态组合，M-1是因为当前符号由状态和输入共同决定。expected_phase预计算每个状态在假设下一个符号为+1时的期望相位，这样在循环里只需要判断bit=0还是1，避免重复计算。

分支度量用的是欧氏距离的平方，sum(abs(rx_seg - exp(1j*exp_phase)).^2)，这是加性高斯白噪声信道下的最优度量。状态转移那行next_state = mod(s-1, 2^(M-2)) * 2 + bit + 1是移位寄存器的逻辑：新状态把旧状态左移一位，低位补当前bit。

回溯部分从最小路径度量的状态出发，反向追踪survivor_path，得到最优比特序列。注意回溯出来的bits是倒序的，实际用的时候要fliplr一下。

提示：Viterbi解调的计算量随状态数指数增长，BT=0.3时M=4或5就够了，再大收益很小但耗时翻倍。如果跑起来太慢，先把M降到3试试。

3. 完整链路搭建：从参数配置到误码率曲线

3.1 主脚本结构与参数配置

把调制和解调串起来的主脚本不需要太长，关键是参数要集中管理。我一般把所有可调参数放在脚本开头，后面直接引用变量名，这样改BT或者SNR的时候不用满篇找数字。

% GMSK完整链路仿真主脚本 clear; close all; clc; % ========== 参数配置 ========== BT = 0.3; % 带宽-码元周期积 Tb = 1e-3; % 码元周期(秒) L = 16; % 过采样倍数 span = 4; % 滤波器截断符号数 M = 4; % Viterbi记忆长度 num_bits = 10000; % 仿真比特数 SNR_dB = 0:2:12; % 信噪比范围 % ========== 生成滤波器 ========== h = gaussian_filter(BT, Tb, L, span); % ========== 误码率循环 ========== ber = zeros(size(SNR_dB)); for idx = 1:length(SNR_dB) % 生成随机比特 data = randi([0 1], 1, num_bits); symbols = 2*data - 1; % 映射到-1/+1 % 调制 [tx, ~] = gmsk_modulate(symbols, h, L); % 加噪 rx = awgn(tx, SNR_dB(idx), 'measured'); % 解调 bits_out = gmsk_demodulate(rx, h, L, M); % 计算误码率 ber(idx) = sum(bits_out ~= data) / num_bits; end % ========== 绘图 ========== semilogy(SNR_dB, ber, 'b-o', 'LineWidth', 1.5); grid on; xlabel('SNR (dB)'); ylabel('BER'); title(sprintf('GMSK误码率曲线 (BT=%.1f, M=%d)', BT, M));

awgn的'measured'选项很重要，它会先测量信号功率再按SNR加噪。如果不加这个选项，MATLAB默认假设信号功率为0dBW，加出来的噪声量级可能完全不对。num_bits取10000是折中，再少误码率曲线会抖，再多跑得慢。

3.2 误码率曲线与理论值对比

跑完主脚本会得到一条BER曲线。BT=0.3、M=4的时候，SNR=8dB左右BER应该降到10^-3量级。如果曲线明显偏离，先检查三个地方：滤波器归一化有没有做、cumsum之前有没有上采样、Viterbi回溯有没有倒序。

和理论MSK的BER对比也有意义。GMSK在BT=0.5时性能和MSK接近，BT=0.3时会有1-2dB的损失，这是频谱效率换来的代价。如果仿真结果比理论值差超过3dB，大概率是Viterbi状态数不够或者滤波器截断太短。

% 理论MSK误码率(相干解调) SNR_lin = 10.^(SNR_dB/10); ber_theory = 0.5 * erfc(sqrt(SNR_lin)); hold on; semilogy(SNR_dB, ber_theory, 'r--', 'LineWidth', 1.5); legend('GMSK仿真', 'MSK理论');

erfc是互补误差函数，MSK相干解调的理论BER就是0.5erfc(sqrt(Eb/N0))。注意这里的SNR是Eb/N0还是Es/N0，awgn默认按Es/N0加噪，如果要用Eb/N0需要减去10log10(L)的修正。

3.3 眼图与相位轨迹观察

误码率曲线是宏观指标，调试的时候还需要看微观波形。眼图和相位轨迹是两个最直观的工具。眼图看码间干扰，相位轨迹看调制是否连续。

% 眼图 eyediagram(real(tx(1:2000)), L*2); title('GMSK I路眼图'); % 相位轨迹 figure; plot(unwrap(angle(tx(1:500)))); xlabel('采样点'); ylabel('相位(rad)'); title('GMSK相位轨迹'); grid on;

眼图张开度好的话，采样时刻的判决门限清晰。如果眼图闭合，说明BT太小导致码间干扰严重，或者滤波器截断不够。相位轨迹应该是分段线性的，每过一个码元周期斜率变化一次。如果轨迹有跳变，检查cumsum之前的数据类型是不是复数。

注意：eyediagram需要Communications Toolbox，如果没有这个工具箱，可以用reshape手动画。把信号按2L长度分段，每段画一条线，叠加起来就是眼图。

4. 避坑与排查：GMSK仿真里最容易翻车的五个地方

4.1 误码率曲线不下降，始终在0.5附近

现象：跑完主脚本，BER曲线是一条水平线，不管SNR怎么增加都不降。

原因：最常见的是符号映射反了。randi([0 1])生成的是0/1，如果直接送进调制器而没有做2*data-1的映射，相位累积会完全错误。另一个可能是Viterbi解调的状态转移矩阵写错了，导致所有路径的度量都差不多。

解决：在调制前打印前10个符号，确认是-1/+1而不是0/1。在Viterbi循环里打印path_metric的最小值，如果所有状态的度量都相等，说明分支度量计算有问题。

4.2 相位轨迹发散或跳变

现象：plot(unwrap(angle(tx)))出来的曲线不是分段线性，而是有突然的跳变或者持续上升。

原因：cumsum之前没有做上采样，或者上采样倍数和滤波器长度不匹配。另一个可能是conv的'same'选项导致输出长度和输入不一致，后续索引越界。

解决：检查upsample(data, L)之后的长度是不是length(data)*L，conv之后是不是用'same'截断到相同长度。如果相位持续上升，说明累积没有归零机制，这是正常的，GMSK相位本来就是无界累积的，只要斜率变化正确就行。

4.3 Viterbi解调速度极慢

现象：num_bits=10000跑了好几分钟还没结束。

原因：状态数太多。M=5时状态数是16，M=6时是32，每增加1状态数翻倍。如果BT=0.5，M=3就够了，没必要用M=5。

解决：先把M降到3，看BER曲线是否可接受。如果必须用大M，把num_bits降到2000先调通，再逐步增加。另外expected_phase的预计算可以移到循环外面，避免重复计算。

4.4 加噪后信号功率不对

现象：awgn加噪后，信号的星座图完全散掉，或者SNR设置和实际不符。

原因：awgn默认假设输入信号功率为0dBW，如果信号幅度不是1，加噪量级会错。'measured'选项可以自动测量，但要求信号是复数且长度足够。

解决：在awgn之前先归一化信号功率：tx = tx / sqrt(mean(abs(tx).^2));。这样信号功率为1，awgn的SNR就是准确的。或者直接用'measured'选项，但注意它只测量输入信号的功率，不测量噪声功率。

4.5 误码率比理论值差很多

现象：BER曲线趋势对，但整体比MSK理论值差3dB以上。

原因：滤波器截断太短导致频谱泄漏，或者Viterbi状态数不够导致路径度量不准。另一个可能是span取太小，高斯滤波器的尾巴被截掉了。

解决：把span从3增加到4或5，看BER有没有改善。如果改善明显，说明截断是瓶颈。如果没改善，把M从3增加到4。两个都试过还不行，检查gaussian_filter里的归一化是不是sum(h)而不是trapz，离散求和和积分的归一化系数不同。

5. 二次开发与参数调优：把这套代码改成你自己的

这套代码的扩展性主要在两个地方：调制端和解调端。调制端可以换滤波器类型，比如把高斯滤波器换成根升余弦，就变成了另一种调制方式。解调端可以换判决算法，比如把Viterbi换成BCJR，得到软输出用于Turbo码迭代译码。

先说参数调优。BT是最关键的参数，它直接决定频谱效率和误码率的权衡。我一般会跑一组BT值：0.2、0.3、0.5、0.7，画在同一张图上对比。BT=0.2时频谱最紧凑但BER最差，BT=0.7时接近MSK性能。如果你的应用对带外泄漏要求严格，选0.3；如果对灵敏度要求高，选0.5以上。

% BT参数扫描 BT_values = [0.2, 0.3, 0.5, 0.7]; figure; for bt = BT_values h = gaussian_filter(bt, Tb, L, span); % ... 跑BER循环 ... semilogy(SNR_dB, ber, 'LineWidth', 1.5); hold on; end legend(arrayfun(@(x) sprintf('BT=%.1f', x), BT_values, 'UniformOutput', false));

再说解调端的改造。Viterbi解调的输出是硬判决比特，如果要接信道译码器，需要软信息。把分支度量改成对数似然比，输出每个比特的LLR，就能和LDPC或Turbo码级联。改动量不大，主要是把branch_metric从欧氏距离改成log(sum(exp(-d^2/2/sigma^2)))的形式。

还有一个实用的改造是加同步。实际接收的信号有频偏和相偏，这套代码假设理想同步。加一个Costas环或者基于前导码的频偏估计，就能处理非理想情况。我一般会在帧头加32个已知符号做前导，接收端用相关峰做定时同步，再用相位差做频偏估计。

最后说一个我踩过的坑：MATLAB的conv函数在信号很长时效率不高，num_bits超过50000会明显变慢。换成filter函数或者用FFT做频域卷积，速度能快3-5倍。如果要做实时仿真，还得把Viterbi循环向量化，用矩阵运算代替逐状态循环。

从那以后我每次改完参数，都会先跑num_bits=1000的快速验证，确认BER趋势对了再上10000的全量仿真。这个习惯帮我省了不少等待时间。希望帮到你。

本文还有配套的精品资源，点击获取