kuopinqpsk.m

本程序对于 循环统计量中存在的直扩信号进行检测

function [jielun]=kuopinQPSK(SNR)
% % %基于循环统计量的直扩信号检测  
% clear all;
% close all;
% clc
% whj

%-----------产生直扩BPSK信号---------------------------
fcarr = 25e3;           % Carrier frequency
fs = 200e3;		        % Sampling frequency
Ts = 1/fs;	            % Sampling time = 1/fs
fm = 2.5e3;               % 信息码速率（波特率） 
Tm = 1/fm;              % Symbol time 码元周期

M_mm=2;%BPSK调制
B_length=100;%符号个数
N_mm=fs/fm;%一个码元周期内含有的采样点数
sigma_A=1;%信号幅度
n_mm=B_length*N_mm;%信号的总长度
symbols_mm=randint(1,B_length,[0,1]);%产生基带码元/
symbolsmm=symbols_mm;
x1_B = ones(N_mm,1)*symbolsmm;	
x1_BB = x1_B(:)'; 
% subplot(511)
% plot(x1_BB)
% axis([0 length(x1_BB) -0.2 3.2])
x1_qpsk=sigma_A*dmod(symbolsmm,fcarr,fm,fs,'psk',M_mm);%BPSK调制信号
% subplot(512)
% plot(x1_qpsk)
% axis([0 length(x1_qpsk) -1.2 1.2])

fs_mm=200e3;%m序列采样率
fm_mm=12.5e3;%波特率  Rc表示
n=length(x1_qpsk);
mm=30;
symbols_rand=idinput(mm,'prbs',[0,1],[-1,1]);
nn=ceil(fs/fm*B_length/mm);
symbols_rand=repmat(symbols_rand,1,nn);
s_rand=symbols_rand(:);
s_rand=repmat(s_rand,1,fs_mm/fm_mm);
s_rand=s_rand';
s_rand=s_rand(:)';%PN码
% subplot(513)
% plot(s_rand)
% axis([0 length(x1_qpsk) -1.2 1.2])

xx_rand=x1_qpsk.*s_rand(1:length(x1_qpsk));  %BPSK调制后扩频
T=1024;
xx_rand=xx_rand(1:T);
% % subplot(514);plot(xx_rand)
% % axis([0 length(x1_qpsk) -1.2 1.2])
% plot(xx_rand)

% %----------高斯白噪声----------------------------------------
% %  SNR=2;     %与正常信号的信噪比
% noise=sqrt(1/(10^(SNR/10)));  %信号幅度取1
% noise=noise*randn(1,length(xx_rand));
% xx_rand=xx_rand+noise;

%-----------alpha低阶稳定分布噪声------------------------------------------
%  SNR =2 ;
% ne1 = randn(1,len);
ne1 = mystable(length(xx_rand),0.9,0)';%设定低阶稳定分布噪声阶数，注意不要和循环频率混了
na = (mse(xx_rand)/10^(SNR/10))*ne1;
xx_rand = xx_rand+na;


%-----------构造检验序列---------------------------
m=40;
df=100;  %间隔步长因子
alpha=[-3*fcarr:fcarr/df:3*fcarr];%循环频率范围
% %%%%%%%%%%%%%%
% xx_rand =(abs(xx_rand).^(-0.5)).*conj(xx_rand);      %低阶化
% %%%%%%%%%%%%%%%
Ruv=csdr(xx_rand,xx_rand,alpha,fs,m);
% plot(abs(Ruv))
Ruv=Ruv';
% gama=0.01;
for k=1:length(alpha)
    Ruv2=[real(Ruv(:,k)),imag(Ruv(:,k))];
   %-----------求协方差矩阵----------- 
   X=Ruv2;
    [n,p]=size(X); 
   X=X-ones(n,1)*mean(X); 
   Y=X'*X/(n-1);
    %------------------------
    Yr=(T*Ruv2*pinv(Y))';        %构造检验统计量(此处可能有错误，可按照文献重新修改运行看结果)
    Prb=0;
    Prs=0;
    Pfa=0.01;                   %设置虚警概率           
   gama=9.210;
%     Pfa=0.001;                   %设置虚警概率           
%     gama=13.815;
   for l=1:length(Yr)
       if(Yr(:,l)>gama)Prb=Prb+1;
       else Prs=Prs+1;
       end
   end
end
   Prb=Prb/length(Yr);
   Prs=Prs/length(Yr);
   if(Prb>Pfa)
       if((k-((length(alpha))-1)/2)~=0)
%            jielun='有直扩信号'
             jielun=1;
%              xunhuanf=(k-((length(alpha))-1)/2)*fcarr/df;%通过估计循环频率来估计载频
       end
   else jielun=0;%jielun='无直扩信号'
   end