无公式搞懂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);
    else
        J_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*n
           if(j<=(i-1)*fn)
               c(j)=0;
           elseif(j<=(i+L-1)*fn)
               c(j)=J_g(j-(i-1)*fn);
           else
               c(j)=0.5;
           end
       end
    else
        for k=1:fn*n
           if(k<=(i-1)*fn)
               c(k)=0;
           elseif(k<=(i+L-1)*fn)
               c(k)=-J_g(k-(i-1)*fn);
           else
               c(k)=-0.5;
           end
       end
    end
    title("每个码元对相位的影响");
    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;
    else
        J_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_g
    Alpha((k-1)*B_sample+j)=pi*J_Alpha+L*pi/2;% 刷新 Alpha
end

k=2;%计算第2个码元的相位
L=0;
for j=1:B_sample%采样点数
    %34码元对应乘以J_g
    J_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;
end

k=3;%计算第3个码元的相位
L=0;
for j=1:B_sample%采样点数
    %345码元对应乘以J_g
    J_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;
end

k=4;%计算第4个码元的相位
L=0;%level为0,表示二进制相位轨迹,因而单位递增pi
for j=1:B_sample%采样点数
    %3456码元对应乘以J_g
    J_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;
end

L=0;%level为0,表示二进制相位轨迹,因而单位递增pi
for k=5:code_len-2
    if k==5
        L=0;
    else
        L=L+code(k-3);  
    end
    for j=1:B_sample%采样点数
        %连续5个码元对应乘以J_g
        J_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     3

function [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;
    else
        modu_signal_shift(i) = modu_signal(i-Dt_c_num/4);
    end
end

xt = 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 = 512
        bt(i)=-1;
    else
        bt(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 Frequency

Fpass = 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]

参考博客:
[1] https://blog.csdn.net/yundanfengqing_nuc/article/details/70211287
[2] https://blog.csdn.net/u013346007/article/details/62917617

  • 8
    点赞
  • 53
    收藏
  • 打赏
    打赏
  • 1
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
©️2022 CSDN 皮肤主题:猿与汪的秘密 设计师:我叫白小胖 返回首页
评论 1

打赏作者

时空默契

你的鼓励将是我创作的最大动力

¥2 ¥4 ¥6 ¥10 ¥20
输入1-500的整数
余额支付 (余额:-- )
扫码支付
扫码支付:¥2
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值