无公式搞懂GMSK调制原理，附详细注释的matlab GMSK调制解调原理仿真源码

一、原理
GMSK调制是一种连续相位调制，每个码元会对信号产生0.5个周期也就是π的相位变化，但单个码元引起的相位变化的过程不在一个码元时间内完成，而是分散到L个码元时间内。假设每个码元对相位的影响在6个码元时间内完成，相位成形函数（高斯曲线的积分）如下图：
在这里插入图片描述
码元对相位的影响按照相位成形函数变化，在6个码元时间内对信号产生 0.5 或 -0.5 个周期的影响。

% 高斯曲线
function J_g = gauss_wave(n)
peak = 1;       %波峰值 peak=1/(sqrt(2*pi)*sigma)
sigma=(1/(peak*sqrt(2*pi)));
a=sigma*6;      %对称轴，离对称轴3*sigma的值就已经接近于0了
%n = 320;        %采样点数=每个码元的采样点数*关联码元的数量
x=0:a*2/n:a*2-a*2/n;
gaussf=peak*exp(-((x-a).^2)/(2*sigma^2));
% figure();
% subplot(211)
% plot(x,gaussf);%高斯曲线积分得到相位成形函数
J_g=zeros(1,length(gaussf));
for i=1:length(gaussf)if i==1 J_g(i)=gaussf(i)*(a*2/n);elseJ_g(i)=J_g(i-1)+gaussf(i)*(a*2/n);end
end
% 
% subplot(212)
% plot(x,J_g);
J_g=J_g./2;

相邻两个码元对信号相位的影响相差一个码元时间，每个码元对相位的影响如下：
在这里插入图片描述
码元为 [1 1 1 0 1 0 1 0 0 0 0 1 1 0]
上图显示了14个码元分别对信号相位影响的时间，将这些码元对相位的影响叠加到一起就得到了下图所示的总的相位变化。将总的相位变化去控制载波的相位就得到了GSMK调制信号了。单个码元对相位的影响在发送前影响为0，在L个码元的时间后为固定的 0.5或 -0.5个周期，对相位的影响只在L个码元时间内变化[1]。
在这里插入图片描述

code = [1 1 1 -1 1 -1 1 -1 -1 -1 -1 1 1 -1];%码元
n = length(code);  %码元数量
fn= 64;         %每个码元采样点数
L = 5;          %码元关联数量
x = 1:fn*n;     %总的采样数
J_g=gauss_wave(fn*L);%相位变化函数
figure(1)
plot(1:1/fn:L+1-1/fn,J_g);
title("相位成形函数");c = zeros(1,fn*n);
p = zeros(1,fn*n);
figure(2)
for i=1:n%第i个码元对相位的影响if(code(i)==1)for j=1:fn*nif(j<=(i-1)*fn)c(j)=0;elseif(j<=(i+L-1)*fn)c(j)=J_g(j-(i-1)*fn);elsec(j)=0.5;endendelsefor k=1:fn*nif(k<=(i-1)*fn)c(k)=0;elseif(k<=(i+L-1)*fn)c(k)=-J_g(k-(i-1)*fn);elsec(k)=-0.5;endendendtitle("每个码元对相位的影响");plot(x/64,c+i-1);
%     if(code(i)==1)
%         axis([0 n -0.1 0.6]);
%     else
%         axis([0 n -0.6 0.1]);
%     end
%     title(['第',num2str(i),'个码元对相位的影响']);p = p + c;hold on;
end
%     subplot(n+1,1,n+1)figure()plot(x/64,p);title("总的相位变化");

二、调制解调matlab仿真源码[2]

调制

%产生GMSK调制信号%输入:非极性码，码元频率，载波频率,采样频率
function modu_signal = b_modulation(code,F_code,F_carry,Fs)%code = [-1 -1 1 1 -1 1 1 -1 -1 1 -1 1 1 -1 1 -1 1 -1 -1 1 -1 1 -1 1 1 -1 -1 1 -1 1 -1 1 -1 -1 1 -1 1 -1 -1 1 1 ];
%code = [-1 -1 1 1 -1 1 1 -1 -1 -1 -1 -1 -1 -1 -1 -1 1 -1 -1 -1 -1 1 -1 1 1 -1 -1 1 1 1 1 1 1 1 1 1 1 -1 -1 1 1 ];
%code = [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1];
%code = [-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1];
%F_code = 16000;
%F_carry = 32000;
%Fs = 32000*4;code_len = length(code);    %码元数量
Tb = 1/F_code;              %码元周期
Dt=1/Fs;                    %采样间隔
B_sample=Tb/Dt;             %每个码元采样点数
code_exp = code_expand(code,B_sample);    %code扩展B_sample倍
t=0:Dt:code_len*Tb-Dt;      %仿真时间离散点
BbTb=0.5;                   %取调制指数为0.5
Bb=BbTb/Tb;                 %3dB带宽
%************************* Filter initialization **************************tt=-2.5*Tb:Dt:2.5*Tb-Dt;   %取5个码元周期进行类卷积gaussf=erfc(2*pi*Bb*(tt-Tb/2)/sqrt(log(2))/sqrt(2))/2-erfc(2*pi*Bb*(tt+Tb/2)/sqrt(log(2))/sqrt(2))/2; J_g=zeros(1,length(gaussf)); %高斯滤波器的积分
for i=1:length(gaussf)if i==1 J_g(i)=gaussf(i)*Dt;elseJ_g(i)=J_g(i-1)+gaussf(i)*Dt;end
end
J_g=J_g/2/Tb; %***************************** 相位计算 ***********************************%计算相位Alpha
Alpha=zeros(1,length(code_exp));k=1;  %计算第1个码元的相位
L=0;
for j=1:B_sample  %每个码元采样点数J_Alpha=code(k+2)*J_g(j);  %第3码元对应乘以J_gAlpha((k-1)*B_sample+j)=pi*J_Alpha+L*pi/2;% 刷新 Alpha
endk=2;%计算第2个码元的相位
L=0;
for j=1:B_sample%采样点数%第3、4码元对应乘以J_gJ_Alpha=code(k+2)*J_g(j)+code(k+1)*J_g(j+B_sample);Alpha((k-1)*B_sample+j)=pi*J_Alpha+L*pi/2;
endk=3;%计算第3个码元的相位
L=0;
for j=1:B_sample%采样点数%第3、4、5码元对应乘以J_gJ_Alpha=code(k+2)*J_g(j)+code(k+1)*J_g(j+B_sample)+code(k)*J_g(j+2*B_sample);Alpha((k-1)*B_sample+j)=pi*J_Alpha+L*pi/2;
endk=4;%计算第4个码元的相位
L=0;%level为0，表示二进制相位轨迹，因而单位递增pi
for j=1:B_sample%采样点数%第3、4、5、6码元对应乘以J_gJ_Alpha=code(k+2)*J_g(j)+code(k+1)*J_g(j+B_sample)+code(k)*J_g(j+2*B_sample)+code(k-1)*J_g(j+3*B_sample);Alpha((k-1)*B_sample+j)=pi*J_Alpha+L*pi/2;
endL=0;%level为0，表示二进制相位轨迹，因而单位递增pi
for k=5:code_len-2if k==5L=0;elseL=L+code(k-3);  endfor j=1:B_sample%采样点数%连续5个码元对应乘以J_gJ_Alpha=code(k+2)*J_g(j)+...code(k+1)*J_g(j+1*B_sample)+...code(k  )*J_g(j+2*B_sample)+...code(k-1)*J_g(j+3*B_sample)+...code(k-2)*J_g(j+4*B_sample);Alpha((k-1)*B_sample+j)=pi*J_Alpha+mod(L,4)*pi/2;end
end%B_num-1;
k=code_len-1;
L=L+code(k-3); %level=0  此时L=1+Ak(4)=1-1=0
for j=1:B_sample%采样点数J_Alpha=code(k+1)*J_g(j+B_sample)+code(k)*J_g(j+2*B_sample)+code(k-1)*J_g(j+3*B_sample)+code(k-2)*J_g(j+4*B_sample);Alpha((k-1)*B_sample+j)=pi*J_Alpha+mod(L,4)*pi/2;%相位变化+0*pi/2 刷新 Alpha512 中 64*6+1-64*7
end%B_num;
k=code_len;%计算最后一个码元的相位
L=L+code(k-3);%level=1  此时L=0+Ak(4)=0+1=1
for j=1:B_sample%采样点数J_Alpha=code(k)*J_g(j+2*B_sample)+code(k-1)*J_g(j+3*B_sample)+code(k-2)*J_g(j+4*B_sample);Alpha((k-1)*B_sample+j)=pi*J_Alpha+mod(L,4)*pi/2;%相位变化+1*pi/2
end  %**************************************************************************
figure()
x=0;
subplot(311)
plot(t/Tb,code_exp,t/Tb,x,'r:');
axis([0 code_len -1.5 1.5]);
title('原始信号');subplot(312)
plot(t/Tb,Alpha*2/pi);
axis([0 code_len min(Alpha*2/pi)-1 max(Alpha*2/pi)+1]);
title('相位波形');S_Gmsk=cos(2*pi*F_carry*t+Alpha);%载波频率+相位
subplot(313)
plot(t/Tb,S_Gmsk,t/Tb,x,'r:');
axis([0 code_len -1.5 1.5]);
title('GMSK波形');modu_signal = S_Gmsk;%频谱
figure();
plot_fq(Fs,modu_signal);
title('modusignal频谱');

码元扩展

% s = code_expand([1 2 3],2)
% s = 1     1     2     2     3     3function [out]=code_expand(d,M)   
N=length(d);           %基带信号码元长度
z=zeros(M,N);          %矩阵M为采样点  N为基带信号码元数量
z(1,:)=d;              %将零矩阵第一行替换成基带信号中的8个码元
zero_replace=reshape(z,1,M*N);  % 1行 m*n 列
temp = conv(zero_replace,ones(1,M));
out = temp(1:N*M);

绘制频谱

%绘制频谱图
function f_img = plot_fq(Fs,signal)
N = length(signal);
n = 1 : N;
f = n.*(Fs/N);
y = abs(fft(signal));
fx = f-f(length(f)/2);%移到原点对称
yy = [y(length(y)/2+1:length(y)) y(1:length(y)/2)];%原点对称
f_img = plot(fx,yy);

解调

%GSMK非相干解调
clear all;
%输入：待解调信号、码元速率、载波频率、采样频率
%function de_siganl = b_demodulation(modu_signal,F_code,F_carry,Fs)
code = [0 0 1 1 -1 1 1 -1 -1 1 -1 1 1 -1 1 -1 1 -1 -1 1 -1 1 -1 1 1 -1 -1 1 -1 1 -1 1 -1 -1 1 -1 1 -1 -1 1 1 ];
%code = [-1 -1 1 1 -1 1 1 -1 -1 -1 -1 -1 -1 -1 -1 -1 1 -1 -1 -1 -1 1 -1 1 1 -1 -1 1 1 1 1 1 1 1 1 1 1 -1 -1 1 1 ];
%code = [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1];
%code = [-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1];
F_code = 1600;
F_carry = 32000;
Fs = 32000*4;%产生调制信号用于解调
modu_signal = b_modulation(code,F_code,F_carry,Fs);Tb = 1/F_code;              %码元周期
Dt=1/Fs;                    %采样间隔
B_sample=Tb/Dt;             %每个码元采样点数
Dt_c=1/F_carry;             %载波周期
Dt_c_num = Dt_c/Dt;         %每个载波周期采样点数%载波相移pi/2
modu_signal_shift = zeros(1,length(modu_signal));
for i=1:length(modu_signal)if(i<=Dt_c_num/4)modu_signal_shift(i) = 0;elsemodu_signal_shift(i) = modu_signal(i-Dt_c_num/4);end
endxt = modu_signal .* modu_signal_shift;%频谱
figure();
plot_fq(Fs,xt);
title("自相干频谱")%低通滤波
Fpass = F_code;
Fstop = 2*F_code;
lp_fir = demodu_lp_fir(Fs,Fpass,Fstop);
y=filter(lp_fir,xt);
figure()
plot(y);
title("经过低通滤波器后波形");for i=1:length(code)if y(i*B_sample)>0  % 8*64 = 512bt(i)=-1;elsebt(i)=1;end
end
figure();
subplot(211)
plot(code);
title("原始信号");
subplot(212);
plot(bt);
title("解调信号");

低通滤波

function Hd = demodu_lp_fir(Fs,Fpass,Fstop)
%DEMODU_LP_FIR Returns a discrete-time filter object.% MATLAB Code
% Generated by MATLAB(R) 9.5 and Signal Processing Toolbox 8.1.
% Generated on: 27-Oct-2021 15:13:41% Equiripple Lowpass filter designed using the FIRPM function.% All frequency values are in MHz.
Fs = Fs/1000;  % Sampling FrequencyFpass = Fpass/1000;              % Passband Frequency
Fstop = Fstop/1000;              % Stopband Frequency
Dpass = 0.057501127785;  % Passband Ripple
Dstop = 0.0001;          % Stopband Attenuation
dens  = 20;              % Density Factor% Calculate the order from the parameters using FIRPMORD.
[N, Fo, Ao, W] = firpmord([Fpass, Fstop]/(Fs/2), [1 0], [Dpass, Dstop]);% Calculate the coefficients using the FIRPM function.
b  = firpm(N, Fo, Ao, W, {dens});
Hd = dfilt.dffir(b);% [EOF]