数字集群MDR中4FSK的matlab仿真

1. 基本概念

物理层的定义h=0.27；T是采样周期，1/4800秒。而D=3h/2T=1.944kHz。D的物理含义是Maximum deviation，最大频偏。

4FSK是恒包络调制，所以在终端功率上，应该比较省电。

4进制的随机信号产生，并把信号强度调整到规范规定的3,1,-1,-3：

M=4;

Ts=1/4800; % 4.8K的波特率；

x=randint(1000,1,M);

x=x.*2-3; %调整到正确的symbol；

plot(x);

部分信号绘制如下：

3210-1-2-3210220230240250260270280290

由于是随机信号，所以是宽带信号。频谱如下。只绘制到信号的采样频率，绘制频谱：

lx=length(x);

Freqs=1/Ts;

freq=[-Freqs/2:Freqs/lx:Freqs/2-Freqs/lx];

S=fftshift(abs(fft(x)));

Sdb=10*log10(S);

plot(freq,Sdb);：

2520151050-5-2500-2000-1500-1000-50005001000150020002500

进行基带滤波，滤波后的信号和频谱为下图。注意信号带宽被限制为根升余弦滤波的3dB带宽以内。这个带宽|cos(f*pi/1920)|=1/2，f>1920，得到f=640+1920=2560Hz。这是一半带宽：

% 滤波器的基本参数：sqrt，0.2。Fs至少要是2Fd，delay用默认的3。

y=rcosflt(x,1,2,'fir/sqrt',0.2,3,0);

%绘制频谱：

ly=length(y);

Freqs=2/Ts;

freq=[-Freqs/2:Freqs/ly:Freqs/2-Freqs/ly];

Sy=fftshift(abs(fft(y)));

Sdby=10*log10(Sy);

figure;

plot(freq,Sdby);

2520152310150-5-10-15-20-5000-4000-3000-2000-10000-1-2-3010002000300040005000

620640660680700720740760

接下去，对4FSK调制进行仿真。

2. 4FSK调制仿真

4FSK有3种方式，方式1，用前面余弦滤波后的基带数据进行频率振荡器的控制。这种方式实现容易，但是频率的稳定性差。方式2，键控法。键控分频器的分频比，以控制2种不同分频脉冲的输出。在经过整形，就可以得到FSK信号。这种在DSP的实现中常用到。方法3，用fskmod函数。

不管用何种方式，DMR规范中，对4FSK的调制的结果给出了要求，要求在给定序列下，输出的调制信号，应该和1.2KHz的sine信号，通过最大偏差2.794KHz的FM调制器输出的信号一致。

方法1： fmmod。

Fmmod是模拟调制函数。它以Fs采样频率来看待输入和输出的数据。DMR系统中，基带信号的带宽小于2*1944=3888Hz。Fc我们可以设计的低一些，以降低采样频率，减少计算量。这样可以设计Fc=15kHz。足够高了。这样，为了让输出的图像细致些，我们设计采样频率Fs=150KHz。这个精度应该购了。

定好了这个调后，回过头看rcosflt函数。这个函数中我们感兴趣的就是它的Fs设计。因为是数字信号，所以rcosflt的Fd可以是1，这个1其实代表了1/Ts=4800Hz。我们的目的是让rcosflt的输出的数据能够和fmmod的Fs=150KHz的采样速率匹配上，于是不难算出，每个输入symbol相当于被采样15k/4800个数据点，也就是31.25个点。取整后，是32个点。反过来，调整采样频率成为32*4800=153.6KHz。

既然采样点数已经确定是32，那么rcosflt的Fs也就确定了，调整为32。所以，matlab代码如下（减少些symbol数目，为20个，便于计算）：

M=4;

Ts=1/4800; % 4.8K的波特率；

x=randint(200,1,M);

x=x.*2-3; %调整到正确的symbol；

% 滤波器的基本参数：sqrt，0.2。Fs至少要是2Fd，取32，delay用默认的3。

y=rcosflt(x,1,32,'fir/sqrt',0.2,3,0); % 过渡补0后，滤波器增益会下降。

y=(3/max(y)).*y; % 调整滤波器增益，恢复信号强度。

% 根升余弦滤波后，输出的就是模拟的基带信号了。

% 模拟FM调制。Fc=15KHz, Fs=153.6KHz, freqdev=3888Hz;（因为滤波后，y最大值就在1左右）

Fc=15000;Fs=153600;freqdev= 648;

ym=fmmod(y,Fc,Fs,freqdev);

% 绘制频谱

lx=length(ym);

Freqs=Fs;

freq=[-Freqs/2:Freqs/lx:Freqs/2-Freqs/lx];

S=fftshift(abs(fft(ym)));

Sdb=10*log10(S);

plot(freq,Sdb);

用上面的方法，根据标准输入，可以看到输出的信号频谱为：

M=4;

Ts=1/4800; % 4.8K的波特率；

for i =1:200

if(mod(i,4)<2) x(i)= 3;

else x(i)=-3;

end

end

% 滤波器的基本参数：sqrt，0.2。Fs至少要是2Fd，取32，delay用默认的3。

y=rcosflt(x,1,32,'fir/sqrt',0.2,3,0); % 过渡补0后，滤波器增益会下降。

y=(3/max(y)).*y; % 调整滤波器增益，恢复信号强度。

% 模拟FM调制。Fc=15KHz, Fs=153.6KHz, freqdev=3888Hz;（因为滤波后，y最大值就在1左右）

Fc=15000;Fs=153600;freqdev= 648;

ym=fmmod(y,Fc,Fs,freqdev);

% 绘制频谱

lx=length(ym);

Freqs=Fs;

freq=[-Freqs/2:Freqs/lx:Freqs/2-Freqs/lx];

S=fftshift(abs(fft(ym)));

Sdb=10*log10(S);

plot(freq,Sdb);

方法2：直接写函数

其实就是模拟FM方法，对输入的信号进行积分。略。

方法3：Fskmod

fskmod(x,M,freq_sep,nsamp)

注意，这个函数是有问题的。我对它进行了修改。

对fskmod.m的代码中修改一行如下：

phaseIncr = (0:nSamp-1)' * (-(M-1):2:(M-1)) * 2*pi * freq_sep/2 * samptime/nSamp;

对fskdemod的代码修改一行如下：

t = [0 : 1/(nSamp*Fs) : 1/Fs - 1/(nSamp*Fs)]';

另外注意，在fskmod这个函数中，已经把x的数据调整到-(M-1):2:(M-1)来进行处理了。

更新后的函数我重新命令为gyzfskmod和gyzfskdemod。另外，这个函数已经把输入的0~M-1调整为-(M-1):2:(M-1)了。符合DMR的需求。

思路还是：

通过gyzfskmod函数得到输出的IQ信号。

对IQ信号进行根升余弦滤波；

进行载波调制。

所以过程如下：

M=4;

Ts=1/4800; % 4.8K的波特率；

% x=randint(100,1,M);

x=[0 1 2 3 2 1 0 1 2 3 2 1 0 ];

Num=length(x);

freqsep=0.27;

nsamp=10;

Fs=1;

% Fc=nsamp*1/Ts; % 载波频率

Fc=10000; % 载波频率

y=gyzfskmod(x,M,freqsep,nsamp,Fs);

注意，这里调制出来的y信号，就是复数表示的模拟IQ信号。对于模拟的IQ信号，不能够直接用数字滤波器rcosflt。正确的如下：

% 设计滤波器，用filterbuilder的Nyquist类型。

% 输出信号的Fs是nsamp/Ts。所以Band=nsamp/2;

% Frequency unit我设置为归一化为1;

% Transition width是 0.2

filterbuilder

%输出Hnyq。这个滤波器的时延没有直接的控制参数。通过解调，可以看到大概时延了2个symbol。

% 输出滤波

yf=filter(Hnyq,y);

% 输入滤波

yfin=filter(Hnyq,y);

% 解调

Z=gyzfskdemod(yfin,M,freqsep,nsamp);

% 输出波形绘图

tlen=length(x)*Ts; % 秒

% 时间分片

t=[0:tlen/length(yf):tlen-tlen/length(yf)];

% 计算载波

ccos=cos(2*pi*Fc.*t);

csin=sin(2*pi*Fc.*t);

% 计算输出波形

S=reshape(real(yf),1,length(yf)).*ccos-reshape(imag(yf),1,length(yf)).*csin;

plot(S);

% 绘制频谱

figure;

p=real(S);

lx=length(p);

Freqs=nsamp/Ts;

freq=[-Freqs/2:Freqs/lx:Freqs/2-Freqs/lx];

St=fftshift(abs(fft(p)));

Sdb=10*log10(St);

plot(freq,Sdb);

1.510.50-0.5-1-1.5020406080100120140