MSK调制解调

1. MSK信号产生原理

1.1 MSK信号时域特征

二进制最小频移键控（Minimum Shift Keying，MSK）信号的标识可写为：
$$S_{MSK}=cos\left({\omega }_{c}t+\frac{\pi a_k}{2T_b}t+{\phi }_k\right)$$
其中，$k{T}_b⩽t⩽(k+1)T_b$，${\omega }_c$是载波角频率，$T_b$是码元宽度，$a_k$是第k个码元中的数据（取值为$\pm 1$）,${\phi }_k$是第k个码元中的相位常数，它在$k{T}_b⩽t⩽(k+1)T_b$中保持不变。
上式中，当$a_k=1$时，信号的频率$f_2=\frac{1}{2\pi }\left({\omega }_c+\frac{\pi }{2T_b}\right)$；当$a_k=-1$时，信号的频率$f_2=\frac{1}{2\pi }\left({\omega }_c-\frac{\pi }{2T_b}\right)$；由此可得频率间隔$△f=f_2-f_1=\frac{1}{2T_b}$，调制指数$h=△f{T}_b=0.5$。
MSK信号和普通的2FSK信号的差别只是选择两个传信频率$f_1$和$f_2$，使这两个频率的信号，在一个码元期间的相位积累严格地相差180°。一般频移键控的两个信号波形具有以下的相关系数：
$$\rho =\frac{sin({\omega }_2-{\omega }_1)T_b}{({\omega }_2-{\omega }_1)T_b}+\frac{sin({\omega }_2+{\omega }_1)T_b}{({\omega }_2+{\omega }_1)T_b}=\frac{sin({\omega }_2-{\omega }_1)T_b}{({\omega }_2-{\omega }_1)T_b}+\frac{sin2{\omega }_c{T}_b}{2{\omega }_c{T}_b}$$
MSK是一种正交调制，其信号波形的相关系数等于零，因此，对MSK信号来说，这个公式等号后面的两项都必须等于零。第一项等于零的条件是$({\omega }_2-{\omega }_1)T_b=2\pi (f_2-f_1)T_b=K\pi (K=1,2,3,\cdots )$，令K等于其最小值1，则$f_2-f_1=\frac{1}{2T_b}$，这正是MSK信号所要求的频率间隔。第二项等于零的条件是$2{\omega }_c{T}_b=4\pi f_c{T}_b=n\pi (n=1,2,3,\cdots )$，即
$$T_b=\frac{n}{4}\frac{1}{f_c}$$
它说明MSK信号在每一个码元周期内，必须包含四分之一载波频率周期的整数倍。
相位常数${\phi }_k$的选择应保证信号相位在码元转换时刻是连续的。根据这一要求，可以导出以下的相位递归条件（也称相位约束条件）。
$${\phi }_k={\phi }_{k-1}+(a_{k-1}-a_k)\frac{\pi k}{2}=\{\begin{array}{cc}{\phi }_{k-1},& a_k=a_k-1\\ {\phi }_{k-1}+k\pi ,& a_k\ne a_k-1\end{array}$$
上式说明，MSK信号在第k个码元的相位常数不仅与当前$a_k$有关，而且与前面的$a_{k-1}$及其相位常数${\phi }_{k-1}$有关。或者说，前后码元之间存在着相关性，对于相干解调来说，${\phi }_k$的起始参考值可以假定为零而不失一般性。因此可以得到${\phi }_k=0或\pi （模2\pi ）$

1.2 MSK信号频谱特性

MSK信号的功率谱密度表示为：
$$W(f)=\frac{16A^2{T}_b}{{\pi }^2}{\left[\frac{cos[2\pi (f-f_c)T_b]}{1-[4(f-f_c)T_b{]}^2}\right]}^2$$
式中，A为信号的振幅。为便于比较，下面直接给出4PSK信号的功率谱密度公式：
$$W(f)=2A^2{T}_b{\left[\frac{sin[2\pi (f-f_c)T_b]}{2\pi (f-f_c)T_b}\right]}^2$$

上图是MSK信号和4PSK信号的功率谱密度曲线。从图中可见，与4PSK信号相比，MSK信号的功率谱密度有较宽的主瓣，前者第一个零点出现于$(f-f_c)T_b=0.5$处，后者第一个零点出现于$(f-f_c)T_b=0.75$处。但是，在主瓣意外，MSK信号的功率谱密度曲线却比4PSK衰减得快。这一点，不难从它们的功率谱密度表示式中看出来，因为当$(f-f_c)T_b\gg 1$以后，4PSK信号的功率谱随$(f-f_c)$的平方而下降，而MSK信号的功率谱却随$(f-f_c)$的四次方而下降。以-3dB电平计算时，MSK信号带宽大于4PSK，以-50dB电平计算和以包含99%的信号功率计算时，MSK信号带宽明显小于4PSK信号。

1.3 MSK信号的产生方法

根据MSK的时域表达式，运用三角函数运算公式，$cos[{\omega }_{c}t+\theta (t)]=cos\theta (t)cos{\omega }_{c}t-sin\theta (t)sin{\omega }_{c}t$，所以MSK信号也可以看成由两个彼此正交的载波$cos{\omega }_{c}t$与$sin{\omega }_{c}t$分别被函数$cos\theta (t)$和$sin\theta (t)$进行振幅调制而合成的。已知$\theta (t)=\frac{\pi a_k}{2T_b}t+{\phi }_k,a_k=\pm 1,{\phi }_k=0或\pi (模2\pi )$，因而
$$cos\theta (t)=cos(\frac{\pi t}{2T_b})cos{\phi }_k,sin\theta (t)=a_{k}tsin(\frac{\pi t}{2T_b})cos{\phi }_k$$
令$cos{\phi }_k=I_k,a_{k}cos{\phi }_k=Q_k$，得到
$$S_{MSK}(t)=I_{k}cos(\frac{\pi t}{2T_b})cos{\omega }_{c}t-Q_{k}sin(\frac{\pi t}{2T_b})sin{\omega }_{c}t$$
通过前面的分析，我们知道，只要找到$I_k、Q_k$与原始数据$a_k$的转换方法，就不难构建产生MSK信号的方法。显然，$I_k、Q_k$的取值只能是1或者-1。
下面给出MSK信号产生的步骤：
(1)对输入数据序列进行差分编码。
(2)把差分编码器的输出数据用串/并转换器分成两路，并且相互交错一个码元宽度$T_b$。
(3)用加权函数$cos(\frac{\pi t}{2T_b})$和$sin(\frac{\pi t}{2T_b})$分别对两路数据进行加权。
(4)用两路加权后的数据分别对正交载波$cos{\omega }_{c}t$和$sin{\omega }_{c}t$进行调制
(5)将两路输出信号进行叠加。

图中左侧虚线由DDS完成MSK数字化调制，具体实现由FPGA完成，右侧为正交调制器AD9977完成正交上变频功能。要求整个调制器输出中频频率为70MHz，如果完全采用FPGA来实现70MHz中频MSK信号的产生，则DDS核必须工作在200MHz以上才能保证调制器输出信号的质量。因此，为了减少FPGA功耗，将MSK分成两个阶段来完成：首先在FPGA中控制DDS完成中心频率为6MHz的MSK调制，输出速率为32MHz的两路正交信号I(t)、Q(t)，然后这两路信号与AD9977里的数字本振64MHz相混频，即完成数字正交调制过程，将中心频率变至70MHz。

2. MSK调制信号的FPGA实现

条件：符号速率Rb=1Mbps，载波信号频率fc=6MHz，采样频率fs=32Rb，输出数据位宽Bout=15。
源码

module MskMod (
	rst,clk,din,
	it,qt);
	
	input    rst;              //复位信号，高电平有效
	input		clk;   				//数据采样时钟:32MHz
	input	   din;              //调制原始数据，Rb=1 MHz
	//输出的CPFSK信号：fc=6 MHz，h=0.5，f1=5.75 MHz，f2=6.25 MHz
	output signed [14:0]	it; 
	output signed [14:0]	qt; 
   
	wire cod;
	//实例化NCO核
	Code	u1 (
		.rst (rst),
		.clk (clk),
		.din (din),
		.dout (cod));

	//实例化频率调制模块核
	FskMod	u2 (
		.rst (rst),
		.clk (clk),
		.din (cod),
		.it (it),
		.qt (qt));

		
endmodule

module Code (
	rst,clk,din,
	dout);
	
	input    rst;              //复位信号，高电平有效
	input		clk;   				//系统时钟信号，与采样频率相同，fs=32 MHz
	input	   din;              //差分编码前的调制数据：1 MHz
	output   dout; 				//差分编码后的调制数据：1 MHz


   reg data;
	reg [4:0] count;
	always @(posedge clk or posedge rst)
		if (rst)
			begin
				data <= 1'b0;
				count <= 5'd0;
			end
		else
			begin
				count <= count + 5'd1;
				if (count==5'd31)
					if (din)
						data <= ~data;
			end
			
	assign dout = data;
		
endmodule

module FskMod (
	rst,clk,din,
	it,qt);
	
	input    rst;              //复位信号，高电平有效
	input		clk;   				//数据采样时钟:32MHz
	input	   din;              //调制原始数据，Rb=1 MHz
	//输出的CPFSK信号：fc=6 MHz，h=0.5，f1=5.75 MHz，f2=6.25 MHz
	output signed [14:0]	it; 
	output signed [14:0]	qt; 


	
	//实例化NCO核所需的接口信号
	wire reset_n,out_valid,clken;
	wire [24:0] carrier;
	wire signed [24:0] frequency_df;
	assign reset_n = !rst;
	assign clken = 1'b1;
	assign carrier=25'd6291456;//6MHz
	
	//实例化NCO核
	dds	u0 (
		.phi_inc_i (carrier),
		.clk (clk),
		.reset_n (reset_n),
		.clken (clken),
		.freq_mod_i (frequency_df),
		.fsin_o (it),
		.fcos_o (qt),
		.out_valid (out_valid));

	//Rb=1MHz，h=0.5，df1=-0.25MHz df2=0.25MHz
	//根据输入数据的高电低平，设置不同的频率偏移量
	assign frequency_df = (din)?25'd262144:-25'd262144;

		
endmodule

3. MSK解调原理

3.1 延迟差分解调

一比特延迟差分解调

假设中频滤波器输出信号为：
$$S_{IF}(t)=R(t)cos[{\omega }_{c}t+\phi (t)]+n(t)$$
式中，R(t)为信号的实时包络，n(t)是噪声或干扰。
中频信号和经过延时$T_b$并移相$\pi /2$的中频信号相乘后，得到
$$S_{mult}(t)=R(t)cos[{\omega }_{c}t+\phi (t)]R(t-T_b)cos[{\omega }_c(t-T_b)+\phi (t-T_b)]$$
经过低通滤波器后，此信号变成
$$y(t)=\frac{1}{2}R(t)R(t-T_b)sin[{\omega }_c{T}_b+\Delta \phi (T_b)]$$
式中，$\Delta \phi (T_b)=\phi (t)-\phi (t-T_b)$
根据MSK信号特征，设${\omega }_c{T}_b=k(2\pi )$（k为整数），则
$$y(t)=\frac{1}{2}R(t)R(t-T_b)sin\Delta \phi (T_b)$$
二比特延迟差分解调

由图可知，经过相乘器与低通滤波器的信号为
$$y(t)=\frac{1}{2}R(t)R(t-2T_b)cos[2{\omega }_c{T}_b+\Delta \phi (2T_b)]$$
式中，$\Delta \phi (2T_b)=\phi (t)-\phi (t-2T_b)$
同样，当$2{\omega }_c{T}_b=k(2\pi )$，有
$$y(t)=\frac{1}{2}R(t)R(t-2T_b)cos[\phi (t)-\phi (t-T_b)]cos[\phi (t-T_b)-\phi (t-2T_b)]-sin[\phi (t)-\phi (t-T_b)]sin[\phi (t-T_b)-\phi (t-2T_b)]$$

3.2 平方环相干解调

由于MSK信号的调制指数h=0.5，经过平方处理后，调制指数实际上变为1。调制指数为1的CPFSK信号，其功率谱中存在着离散分量，即二倍传信号频率$2f_H$和二倍空号频率$2f_L$。因此，可以用两个锁相环电路分别提取出这两个频率。因此可以得到载波频率$f_c$和时钟频率$f_R$之间的关系为：
$$f_c=\frac{f_H+f_L}{2},f_R=2f_H-2f_L$$
为了得到$f_c$和$f_R$，在电路中将两个锁相环锁定的$2f_H$和$2f_L$两个信号相乘，并经过低通滤波，就可以得到时钟频率$f_R$信号，它经脉冲形成后得到所需的时钟脉冲，再由此产生各种定时信号；将除以2后得到$f_H$和$f_L$相乘，并讲过低通滤波可获得$f_R/2$，其速率与正交两路数据的速率相同，可用于正交两路数据的判决定时。$2f_H$和$2f_L$除以2后得到的信号为：
$$S_1(t)=cos(2\pi f_{H}t)=cos\left(2\pi f_{c}t+\frac{\pi t}{2T_b}\right)$$
$$S_2(t)=cos(2\pi f_{L}t)=cos\left(2\pi f_{c}t-\frac{\pi t}{2T_b}\right)$$
将上式相加/相减得到
$$S_1(t)+S_2(t)=2cos\frac{\pi t}{2T_b}cos(2\pi f_{c}t)$$
$$S_1(t)-S_2(t)=-2sin\frac{\pi t}{2T_b}sin(2\pi f_{c}t)$$
为便于对比，我们将MSK信号的时域表达式重写如下
$$S_{MSK}(t)=I_{k}cos\left(\frac{\pi t}{2T_b}\right)cos{\omega }_{c}t-Q_{k}sin\left(\frac{\pi t}{2T_b}\right)sin{\omega }_{c}t$$
容易看出，将相干载波信号分别与接收到的中频信号相乘，并滤除相干载波的2倍频信号，经过正确的抽样判决后，即可获得两路正交的调制数据。

4. MSK解调的MATLAB仿真

条件：符号速率Rb=1Mbps，载波信号频率fc=3MHz，采样频率fs=16Rb。
由于锁相环的设计比较复杂，本程序不适用锁相环电路仿真，直接用产生调制信号的载波进行相干解调。
源码

ps=1*10^6;  %码速率为1MHz
Fs=16*10^6; %采样速率为16MHz
fc=3*10^6;  %载波频率为3MHz

N=100;        %数据码元个数
Len=N*Fs/ps;  %仿真数据的长度
x = randint(N,1,2)'; % 产生随机数据做为数据码元
%x=ones(1,N);
dx=ones(1,N);
for i=1:N
    if x(i)==0
        x(i)=-1;
    end
end
%求原码的相对码dx
for i=2:N
    if x(i)==1
        dx(i)=-dx(i-1);
    else
        dx(i)=dx(i-1);
    end
end

%将相对码按奇偶序号分成两路数据，形成Ik\Qk
di=ones(1,N);
dq=ones(1,N);
%取dx的偶数位，并列两位为di
for i=2:2:N
    di(i:i+1)=dx(i);
end
%取dx的奇数位，并列两位为dq
for i=1:2:N-1
    dq(i:i+1)=dx(i);
end

%对原始BIT数据进行Fs/ps倍重采样
udi=ones(1,N*Fs/ps);
udq=ones(1,N*Fs/ps);
for i=1:N
    udi(Fs/ps*(i-1)+1:Fs/ps*i)=di(i);
    udq(Fs/ps*(i-1)+1:Fs/ps*i)=dq(i);
end
% udi=rectpulse(di,Fs/ps);
% udq=rectpulse(dq,Fs/ps);

%产生MSK信号所需的载波信号
t=0:1/Fs:(Len-1)/Fs;
cf0c=cos(2*pi*fc.*t);
sf0c=sin(2*pi*fc.*t);
cfps=cos(pi*ps/2.*t);
sfps=sin(pi*ps/2.*t);

%正交调制法产生msk信号
msk=udi.*cfps.*cf0c-udq.*sfps.*sf0c;

 
%MSK 解调
%第一种相干解调，需要同时获取fc,fb的载波频率
% demod_i=msk.*cf0c.*cfps;
% demod_q=msk.*sf0c.*sfps;
%平方环相干解调，只需获取fL、fH的载波频率
fL=cos(2*pi*fc.*t-2*pi*ps/4.*t);%cos(2*pi*ps/4.*t)
fH=cos(2*pi*fc.*t+2*pi*ps/4.*t);

demod_i=msk.*(fH+fL);
demod_q=msk.*(fH-fL);
%低通滤波后，获取I、Q支路基带波形
b=fir1(30,0.5*ps*2/Fs);%设计低通滤波器
f_i=filter(b,1,demod_i);
f_q=filter(b,1,demod_q);

%求Tb/2的时钟信号波形
% fL2=cos(4*pi*fc.*t-4*pi*ps/4.*t);%cos(2*pi*ps/4.*t)
% fH2=cos(4*pi*fc.*t+4*pi*ps/4.*t);
% b2=fir1(30,1*ps*2/Fs);%设计低通滤波器
% fb=fL2.*fH2;
% rb=filter(b2,1,fb);

fb=fL.*fH;
% tb=filter(b,1,fb);
rb=filter(b,1,fb);
% 
% for i=1:length(rb)-1
%     rb(i)=tb(i+1)-tb(i);
% end
% rb=rb>=0;
%绘制解后的I、Q支路基带波形及时钟信号波形
figure(1); 
num=800;
x_t=0:1/Fs:(num-1)/Fs;
x_t=x_t*10^6;
subplot(411);plot(x_t,udi(1:num));axis([0 25 -1.2 1.2]);
legend('I支路原始数据时域波形');
xlabel('时间(us)');ylabel('幅度(V)');
subplot(412);plot(x_t,f_i(1:num),'-',x_t,rb(1:num),'--');
legend('I支路解调数据','时钟信号');
xlabel('时间(us)');ylabel('幅度(V)');
subplot(413);plot(x_t,udq(1:num));axis([0 25 -1.2 1.2]);
legend('Q支路原始数据时域波形');
xlabel('时间(us)');ylabel('幅度(V)');
subplot(414);plot(x_t,f_q(1:num),'-',x_t,rb(1:num),'--');
legend('Q支路解调数据','时钟信号');
xlabel('时间(us)');ylabel('幅度(V)');




%绘制MSK信号的频谱及时域波形
figure(2);
%求MSK信号的FFT变换
m_msk=20*log10(abs(fft(msk,2048)));
m_msk=m_msk-max(m_msk);
%设置幅频响应的横坐标单位为MHz
x_f=1:length(m_msk);
x_f=x_f*Fs/length(m_msk)/10^6;
%绘制MSK信号的频谱
subplot(211);
plot(x_f,m_msk);axis([0 Fs/2/10^6 -80 0]);
legend('MSK信号频谱');
xlabel('频率(MHz)');ylabel('幅度(dB)');grid on;

%绘制MSK信号的时域波形
subplot(212);
%设置横坐标单位为us
LT=100;
x_t=0:1/Fs:(LT-1)/Fs;
x_t=x_t*10^6;
plot(x_t,msk(100:100+LT-1))
legend('MSK信号时域波形');
xlabel('时间(us)');ylabel('幅度(V)');

运行结果图如下。从图中可以看出，对于I支路数据而言，1/2倍码速率时钟信号与解调后的I路数据完全同步，最佳判决时刻为时钟波形的波峰出；对于Q支路数据而言，1/2倍码速率时钟信号与解调后的I路数据完全正交，最佳判决时刻为时钟波形的波谷出。