[fh-gfsk] fh-gfsk signal analysis and blind demodulation research

1. Software version

matlab2021a

2. System principle

      In wireless communication , In order to combat channel fading and ensure the safety of transmission information, frequency hopping is often used （Frequency Hopping, FH） Communications technology . Gaussian frequency shift keying （Gauss Frequency Shift Keying, GFSK） With constant amplitude envelope 、 Power spectrum concentration 、 Narrow spectrum and other desirable characteristics of wireless communication systems .FH-GFSK It combines the advantages of the above two technologies , therefore , It is widely used in digital communication . This topic is based on GNU Radio First of all, the software radio receiver of FH-GFSK Signal acquisition , Then analyze the collected signals , Finally, the blind demodulation of the signal is realized . Through the research of this subject , It aims to deepen students' understanding of frequency hopping communication 、 Understanding of Gaussian frequency shift keying technology , Master the analysis and processing methods of communication signals

1. complete 5.8GHz Frequency band FH-GFSK Signal acquisition 、 Analysis and blind demodulation .

2. The adoption is based on GNU Radio Software radio receiver for RF FH-GFSK To collect .

3. use Matlab Analyze the speed hopping and frequency hopping patterns of the collected signals , Get its frequency hopping pattern .

4. use Matlab Yes GFSK Gaussian filtering parameters of the signal 、 The modulation index and other parameters are estimated , Used to complete FH-GFSK Blind demodulation of signal .

According to the references ,GMSK The demodulation structure is as follows ：

May refer to ：https://wenku.baidu.com/view/1631562eb307e87101f69679.html

So when making blind estimation , Need to know fc, And parameters of low-pass filter , Behind it is the fixed structure , There is no need to estimate .

Then in the title ：

Gaussian filter parameters , In fact, it is to obtain the corresponding low-pass filter ;

Then the formula of modulation index is ：

       therefore , Estimate modulation index , It's actually getting fd and Rb, and Rb Indicates the symbol rate , As long as the signal is received , You can get it , Don't estimate , So it is estimated fd.FSK Modulation is a modulation method in which the frequency of the carrier is proportional to the information symbol , When sending a message symbol 1 The transmission frequency moves upward fdHz, When sending a message symbol -1 The transmission frequency moves down fdHz.  Then here is FM communication , So estimate fd. therefore , here , Parameter estimation of busy demodulation , The essence is frequency estimation and filter parameter estimation .

3. Part of the source code

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(' Time s')
ylabel(' frequency 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

% Lock the frequency point 
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
 
% Frequency hopping period 
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));
% Modulation index 
% Gaussian filter parameter estimation 
indx = 0;
for i=0:frameNumber-100% Real time parameter estimation by frame 
    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));
        % Remove the directly combined signal after interruption 
        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);% No signal area 
     end
end
axes(handles.axes2);
[R,C] = size(selectRxFrame2);
Rx    = reshape(selectRxFrame2',[1,R*C]);
plot(real(Rx),'b');
xlabel(' Time s')
grid on;
axis([2000,20000,-1,1]);

% The data is restored to the original one-dimensional signal 
[R,C] = size(selectRxFrame2);
Rx    = reshape(selectRxFrame2',[1,R*C]);

 

% Filter estimate 
% Computational bandwidth , By calculation -3db Spectrum range as bandwidth 
[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(' Modulation index ');
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);
    % Filtering is time-consuming , I only intercept a part here for processing 
    Isignal  = Isignal_(1:length(RR));
    Qsignal  = Qsignal_(1:length(RR));
    % Output binary data 
    tmps = Isignal.*[diff(Qsignal,1),0] - Qsignal.*[diff(Isignal,1),0];
    % Do filtering 
    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(' Skip period ');
cycle
disp(' Frequency point ');
detectFHResult2

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

4. Simulation conclusion

 A01-152