【FH-GFSK】FH-GFSK信号分析与盲解调研究

1.软件版本

matlab2021a

2.系统原理

      在无线通信中，为了对抗信道衰落和保证传输信息的安全常采用跳频（Frequency Hopping, FH）通信技术。高斯移频键控（Gauss Frequency Shift Keying, GFSK）具有恒幅包络、功率谱集中、频谱较窄等无线通信系统所希望的特性。FH-GFSK则结合了上述两种技术的优点，因此，在数字通信中得到广泛应用。本课题采用基于GNU Radio的软件无线电接收机首先对FH-GFSK信号进行采集，随后对采集到的信号进行分析，最终实现信号的盲解调。通过本课题的研究，旨在加深学生对跳频通信、高斯移频键控技术的理解，掌握通信信号的分析与处理方法

1.完成5.8GHz频段的FH-GFSK信号的采集、分析和盲解调。

2.采用基于GNU Radio的软件无线电接收机对射频FH-GFSK进行采集。

3.采用Matlab对采集的信号进行跳速和跳频图案分析，得到其跳频图案。

4.采用Matlab对GFSK信号的高斯滤波参数、调制指数等参数进行估计，用于完成FH-GFSK信号的盲解调。

根据参考文献可知，GMSK解调结构如下所示：

可参考：https://wenku.baidu.com/view/1631562eb307e87101f69679.html

所以在进行盲估计时候，需要知道fc，和低通滤波器的参数，后面的是固定结构，就不用估计了。

然后题目中讲到：

高斯滤波器参数，实际上就是获得对应的低通滤波器；

然后调制指数的计算公式为：

       因此， 估计调制指数，实际上就是获得fd和Rb，而Rb表示符号速率，这个只要信号接收到，就可以得到的，不用估计，所以就是估计fd。FSK调制是载波的频率随信息符号成正比的一种调制方式，当发送信息符号1时发射频率向上搬移fdHz，当发送信息符号-1时发射频率向下搬移fdHz。 然后我们这里是调频通信，因此估计fd。所以，这里，忙解调的参数估计，本质就是就是频率估计和滤波器参数估计两个部分。

3.部分源码

fs          = 100e6;
Nfft        = 4096;
frameNumber = floor(length(x)/Nfft);
txBlockFFT  = zeros(frameNumber,Nfft);
for i = 0:frameNumber-1
    i
    start                   = i*Nfft;
    txBlockFFT(i+1,:)       = fftshift(fft(x(start+1:start+Nfft)));
   [maxValue maxIndex(i+1)] = max(abs(txBlockFFT(i+1,:)));
end
fc          = 5.8e9;
delta_f     = fs/Nfft;
f           = delta_f: delta_f: fs;
f           = f - fs/2;

%detect vaalid signal
validIndexCount             = 1;
validIndex                  = zeros(1,1);
validIndex(validIndexCount) = 0;
validFrameCount = 0;

axes(handles.axes1);
for i = 1:frameNumber-1
    i
    if(max(abs(txBlockFFT(i+1,:))> 280))
        validFrameCount = validFrameCount + 1;
        validFrame(validFrameCount) = i+1;
        if(abs(maxIndex(i+1) - validIndex(validIndexCount)) > 40)
           validIndexCount             = validIndexCount +1;
           validIndex(validIndexCount) = maxIndex(i+1) ;
           detectFHFc                  = validIndex(2:end)*fs/Nfft- fs/2;
        end
        detectFHResult(i+1) = maxIndex(i+1)*fs/Nfft- fs/2+fc;
        t = (0:length(detectFHResult)-1)*Nfft/fs;
        plot(t,detectFHResult,'c*');
        hold on
    end
    axis([0,0.1983,5.74e9,5.86e9]);
    pause(0.001);
    
end

 
hold on;
y_label = fc*ones(1,length(detectFHResult));
plot(t,y_label,'r');
ylim([fc-50e6 fc+50e6]);
xlabel('时间s')
ylabel('频率Hz');
grid on;

ind1 = find(abs(detectFHResult)>0) ;
ind2 = find(detectFHResult==0) ;
detectFHResult(ind2)=[];



flag = [];
for i = 1:length(ind1)-1
    if ind1(i+1)-ind1(i) > 20
       flag = [flag,i]; 
    end
end

%锁定频率点
for i = 1:length(flag)
    if i == 1
       detectFHResult2(i) = mean(detectFHResult(1:flag(i))); 
    else
       detectFHResult2(i) = mean(detectFHResult(flag(i-1)+1:flag(i)));   
    end
end
 
%跳频周期
for i = 1:length(flag)
    if i == 1
       ind12(i) = ind1(flag(i)+1)-ind1(1); 
    else
       ind12(i) = ind1(flag(i)+1)-ind1(flag(i-1)+1); 
    end
end
cycle = floor(mean(ind12));

parameters;
t = (0:Nfft-1)/fs;
t = t';
for i=0:frameNumber-1
    start             = i*Nfft;
    txBlockFFT(i+1,:) = fftshift(fft(x(start+1:start+Nfft)));    
end
Avgs = 1000*mean2(abs(txBlockFFT));
%调制指数
%高斯滤波参数估计
indx = 0;
for i=0:frameNumber-100%分帧进行实时参数估计
    i
    [maxValue,maxIndex(i+1)]= max(abs(txBlockFFT(i+1,:)));
     if max(abs(txBlockFFT(i+1,:)))> Avgs
        indx = indx + 1;
        if(abs(maxIndex(i+1) - validIndex(validIndexCount)) > 40)
           validIndexCount             = validIndexCount +1;
           validIndex(validIndexCount) = maxIndex(i+1) ;
           detectFHFc                  = validIndex(2:end)*fs/Nfft- fs/2;
        end
        detectFHResult(i+1) = maxIndex(i+1)*fs/Nfft- fs/2+fc;
        selectFrameData     = txBlockFFT(i+1,:);
        [maxValue,maxIndex] = max(abs(selectFrameData));
        selectFHFc          = maxIndex*fs/Nfft- fs/2;
        startIndex          = (i+1)*Nfft;
        selectRxFrame(i+1,:) = x(startIndex+1:startIndex+Nfft).*(cos(2*pi*selectFHFc*t)-sqrt(-1)*sin(2*pi*selectFHFc*t));
        %去除中断后的直接合并的信号
        selectRxFrame2(indx,:)= x(startIndex+1:startIndex+Nfft).*(cos(2*pi*selectFHFc*t)-sqrt(-1)*sin(2*pi*selectFHFc*t));
     else
        startIndex          = (i+1)*Nfft; 
        selectRxFrame(i+1,:)= x(startIndex+1:startIndex+Nfft);%无信号区域
     end
end
axes(handles.axes2);
[R,C] = size(selectRxFrame2);
Rx    = reshape(selectRxFrame2',[1,R*C]);
plot(real(Rx),'b');
xlabel('时间s')
grid on;
axis([2000,20000,-1,1]);

%数据还原成原始的一维信号
[R,C] = size(selectRxFrame2);
Rx    = reshape(selectRxFrame2',[1,R*C]);

 

%滤波器估计值
%计算带宽,通过计算-3db频谱范围作为带宽
[x0,t,ssf,yy] = plotspec(Rx,1/fs);
Y2            = 10*log10(yy/max(yy));
Y3            = Y2(length(Y2)/2:end);
indx          = find(Y3>=-3);
ssf2          = ssf(length(ssf)/2:end);
BB            = (ssf2(indx(end))-ssf2(indx(1)));
BT            = 100*BB/fs;
disp('调制指数');
set(handles.edit2,'string',num2str(BT));
Rfinal        = [];
for i = 1:R
    i
    RR           = selectRxFrame2(i,:);
    [Isignal_,h] = glpfsignal(real(RR),fs,BT);
    [Qsignal_,h] = glpfsignal(imag(RR),fs,BT);
    %滤波这块比较费时，我这里只截取一部分进行处理
    Isignal  = Isignal_(1:length(RR));
    Qsignal  = Qsignal_(1:length(RR));
    %输出二进制数据
    tmps = Isignal.*[diff(Qsignal,1),0] - Qsignal.*[diff(Isignal,1),0];
    %做下滤波处理
    tmps = tmps-mean(tmps);
    w    = hamming(128);
    tmps = conv(tmps,w);
    tmps = tmps-mean(tmps);
    Rfinal = [Rfinal,tmps];
end

tmps2= Rfinal>=0;

axes(handles.axes3);
plot(Rfinal)
axis([2000,20000,-2,3]);
axes(handles.axes4);
plot(tmps2)
axis([2000,20000,-1,2]);
 
disp('跳周期');
cycle
disp('频率点');
detectFHResult2

fid = fopen('a.txt','wt');
for i = 1:16
    fprintf(fid,'%6.2f\n',detectFHResult2(i)); 
end
fclose(fid);

4.仿真结论

 A01-152