m分别通过matlab和FPGA实现基于高阶循环谱的信号载波调制识别(四阶循环累量)仿真

1.算法描述

   通信信号调制方式自动识别在军事和民用领域的信号检测、威胁干扰分析、频谱监测管理等方面有着广泛应用,是非合作通信的关键技术。其主要任务是对截获信号进行处理,比如变频、去载波、参数估计等,从中提取特征参数,然后利用适当的算法对调制信号的调制方式做出判决,为后续信号解调等处理提供依据。  

   由于数字调制信号越来越多地应用于通信信号处理领域，因此对数字信号调制识别的研究也越来越多。传统的调制识别的判决方法有：决策判决法、高阶累积量算法和人工神经网络算法等。但是决策判决法在低信噪比环境中识别率不高，而基于人工神经网络的识别方法计算复杂度较高。信号的高阶累积量算法具有很好的抗噪性能，故对基于高阶累积量的通信信号调制识别算法的研究受到了广泛重视。文献利用高阶累积量实现了对 2ASK／BPSK，4ASK，4PSK，2FSK，4FSK信号的分类。文献利用四阶和六阶累积量实现了对 2ASK，4ASK，8ASK，QPSK，8P-SK，16QAM信号的分类。文献利用二、四、六阶累积量实现了对 2ASK／BPSK，4ASK，QPSK，2FSK，4FSK，8FSK，16QAM信号的分类。文献对高阶累积量的四阶、五阶累积量进行了优化和仿真，但是在低信噪比的环境下，对信号的识别率都不高。

   所谓高阶统计量，通常是指信号的高阶矩，高阶累积量等信号统计信息的统称。在这些高阶统计量中，高阶累积量具有十分重要的重用，不同阶数的累积量可以反映出信号的不同特征，因此高阶累积量往往应用在信号分类，信号调制方式识别等领域。高阶累积量中的二阶累积量和四阶累积量可以有效抑制高斯白噪声的干扰，且对相位偏移具有一定的容错能力，其数学表达式为：

下面摘取调制方式在高阶累积量上的差异，

2.仿真效果预览
Quartus II 12.1 (64-Bit)，

ModelSim-Altera 6.6d (Quartus II 11.0) Starter Edition

3.verilog核心程序
matlab部分：

close all;
clear all;
addpath 'func\'
rng('default')
SEL     = 3;
%假定接收端已经实现载波同步
seldata = randi([0 1],1,1000);
fs      = 10000;
R       = 250;
fc      = 1900;
if SEL == 1
   y_2psk1 = func_BPSK(seldata,fs,fc,R,1);
   y_2psk2 = func_BPSK(seldata,fs,fc,R,1.05); 
   com     = y_2psk1+y_2psk2;
end
if SEL == 2
   y_4psk1 = func_QPSK(seldata,fs,fc,R,1);
   y_4psk2 = func_QPSK(seldata,fs,fc,R,1.05);     
   com     = y_4psk1+y_4psk2;
end
if SEL == 3
   y_8psk1 = func_8PSK(seldata,fs,fc,R,1);
   y_8psk2 = func_8PSK(seldata,fs,fc,R,1.05);       
   com     = y_8psk1+y_8psk2;
end
 
%%%%%%%%四阶循环累积量%%%%%%%%%%%%
NOI = max(abs(com))/100;
y   = NOI + awgn(com,30,'measured');    
yft = abs(fft(y,4096)); 
figure;
plot(fs*(0:2047)/4096,yft(1:2048)) 
title('频谱');
fc    = fc;   
alpha =(fc-R):1:(fc+R);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
S4 = y.^4; 
S2 = y.^2; 
m  = size(y); 
k  =(1:m(2))'; 
n  = 0; 
for i=1:(2*R+1)
    i
    e1(:,i)  = exp(-j*8*pi*alpha(i)*k./fs); 
    su1(:,i) = S4'.*e1(:,i); 
    C1(i)    = mean(su1(:,i)); 
    e2(:,i)  = exp(-j*4*pi*alpha(i)*k./fs); 
    su2(:,i) = S2'.*e2(:,i); 
    C2(i)    = mean(su2(:,i)); 
    n        = n+1; 
end 
C=C1-3*C2.^2; 
figure
plot(alpha,abs(C)); 
if SEL == 1
   title('BPSK');
end
if SEL == 2
   title('QPSK');
end
if SEL == 3
   title('8PSK');
end
C40 = C;
 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
S4  = (y.^3).*conj(y); 
S21 = y.^2; 
S22 = y.*conj(y);
m  = size(y); 
k  =(1:m(2))'; 
n  = 0; 
for i=1:(2*R+1)
    i
    e1(:,i)  = exp(-j*8*pi*alpha(i)*k./fs); 
    su1(:,i) = S4'.*e1(:,i); 
    C1(i)    = mean(su1(:,i)); 
    
    e2(:,i)  = exp(-j*4*pi*alpha(i)*k./fs); 
    su2(:,i) = (S21'.*e2(:,i)); 
    C2(i)    = mean(su2(:,i));
    
    e3(:,i)  = exp(-j*4*pi*alpha(i)*k./fs); 
    su3(:,i) = (S22'.*e3(:,i)); 
    C3(i)    = mean(su3(:,i));    
    
    n        = n+1; 
end 
C=C1-3*C2.*C3; 
figure
plot(alpha,abs(C)); 
if SEL == 1
   title('BPSK');
end
if SEL == 2
   title('QPSK');
end
if SEL == 3
   title('8PSK');
end
C41 = C;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
S4  = (y.^2).*(conj(y).^2); 
S21 = y.^2; 
S22 = conj(y).^2;
S3  = y.*conj(y);
m   = size(y); 
k   =(1:m(2))'; 
n   = 0; 
for i=1:(2*R+1)
    i
    e1(:,i)  = exp(-j*8*pi*alpha(i)*k./fs); 
    su1(:,i) = S4'.*e1(:,i); 
    C1(i)    = mean(su1(:,i)); 
    
    e2(:,i)  = exp(-j*4*pi*alpha(i)*k./fs); 
    su2(:,i) = (S3'.*e2(:,i)); 
    C2(i)    = mean(su2(:,i));
       
    e4(:,i)  = exp(-j*4*pi*alpha(i)*k./fs); 
    su4(:,i) = (S22'.*e4(:,i)).*(S22'.*e3(:,i)); 
    C4(i)    = mean(su4(:,i));      
    n        = n+1; 
end 
C=C1 - 2*C2.^2 - C4; 
figure
plot(alpha,abs(C)); 
if SEL == 1
   title('BPSK');
end
if SEL == 2
   title('QPSK');
end
if SEL == 3
   title('8PSK');
end
C42 = C;
FM = round([abs(C41/C40);abs(C41/C42)]);
FM

module C4012_cal(
               i_clk,
                    i_rst,
                    i_din,
                    o_C40,
                    o_FM
              );
                  
input i_clk;
input i_rst;
input signed[15:0]i_din;
output signed[15:0]o_C40;                  
output[1:0]   o_FM;
                  
reg signed[31:0]S2s;
reg signed[63:0]S4s;    
always @(posedge i_clk or posedge i_rst)
begin
     if(i_rst)
      begin
      S2s <= 32'd0;
      S4s <= 64'd0;
      end
else begin
      S2s <= i_din * i_din;
      S4s <= i_din * i_din * i_din * i_din;
     end
end
 
 
wire signed[15:0]S2;
wire signed[15:0]S4;
 
assign S2 = S2s[31:16];
assign S4 = S4s[63:48];
 
 
reg[11:0]cnt;
always @(posedge i_clk or posedge i_rst)
begin
     if(i_rst)
      begin
      cnt <= 12'd0;
      end
else begin
          if(cnt == 12'd2000)
             cnt <= 12'd2000;
      else
           cnt <= cnt + 12'd1;  
     end
end
 
reg start;
always @(posedge i_clk or posedge i_rst)
begin
     if(i_rst)
      begin
      start <= 1'd0;
      end
else begin
          if(cnt == 12'd406)
             start <= 1'd1;
      else
           start <= 1'd0;  
     end
end    
    
    
wire signed[15:0]DOUT;    
C4012 C4012_u( 
              .CLK  (i_clk),
                  .RST  (i_rst),
                  .START(start),
                  .DIN  (i_din[15:6]),
                  .DOUT (DOUT),
                  .FM   (o_FM)
                  );    
                  
 
 
 
reg       [15:0] CNT;
reg signed[15:0]o_C40;
reg             flag;
always @(posedge i_clk or posedge i_rst)
begin
     if(i_rst)
      begin
      flag <= 1'b0;
      end
else begin    
          if(DOUT >= 1000)
             flag <= 1'b1;
      else
          flag <= flag;      
     end
end
 
always @(posedge i_clk or posedge i_rst)
begin
     if(i_rst)
      begin
      CNT <= 16'd0;
      end
else begin    
          if(flag == 1'b1)
             begin
                 CNT <= CNT + 16'd1;
                 if(CNT >= 16'd10000)
                 CNT <= 16'd10000;
             end
     end
end
 
always @(posedge i_clk or posedge i_rst)
begin
     if(i_rst)
      begin
      o_C40 <= 16'd0;
      end
else begin    
         if(CNT <= 100)
           o_C40 <= DOUT;
      else
          o_C40 <= 16'd0; 
     end
end
    
endmodule                   
 01_107m