untitled.m

在工程里识别不同FSK和PSK的三阶循环特性的程序

clear;
clc;
fs=5000    ;         %采样频率
nd=96       ;        %number of data没什么影响
srr=200        ;     %symbol rate
sr=100       ;       %码速率
fc=2000;             %载频1
fc1=2100    ;        %载频2
fc2=2000    ;        %载频3
fc3=2200;            %载频4
ndd=192       ;      %number of data
ss=0;
jielun=0;
 
for xuanhuan=1:1
jieguo=0;

% ********************* signal1********************************
A=0;
B=0;
fff=200;
N=ndd*fs/srr;
xxx=randint(nd,1,4);
M=2;
for i_1=1:nd
    if xxx(i_1)==0
    A=2;
    B=1/2;
elseif xxx(i_1)==1
    A=2;
    B=3/2;
elseif xxx(i_1)==2
    A=2;
    B=-1/2;
elseif xxx(i_1)==3
     A=2;
     B=-3/2;
end
    s_1(i_1)=A*exp(j*2*pi*B*0);                      %产生基带调制相位
    for ii_1=1:fs/sr
        s_2((i_1-1)*fs/sr+ii_1)=s_1(i_1)*exp(j*2*pi*((i_1-1)*fs/sr+ii_1)*(fc+B*fff)/fs);     %产生调制信号
    end
end
% ********************* signal2 ********************************
C=0;
D=0;
N=ndd*fs/srr;
xxx=randint(ndd,1,4);
M=2;
for i_2=1:ndd
    if xxx(i_2)==0
    C=2;
    D=0.5;
elseif xxx(i_2)==1
    C=2;
    D=1;
elseif xxx(i_2)==2
    C=2;
    D=1.5;
elseif xxx(i_2)==3
    C=2;
    D=2;
end
    s_3(i_2)=C*exp(j*2*pi*D/M);
    for ii_2=1:fs/srr
        s_4((i_2-1)*fs/srr+ii_2)=s_3(i_2)*exp(j*2*pi*((i_2-1)*fs/srr+ii_2)*fc1/fs);
    end
end
% ********************* signal3 ********************************
E=0;
F=0;
N=ndd*fs/srr;
xxx=randint(ndd,1,4);
M=2;
for i_3=1:ndd
    if xxx(i_3)==0
    E=2;
    F=1;
elseif xxx(i_3)==1
    E=2;
    F=1;
elseif xxx(i_3)==2
    E=2;
    F=2;
elseif xxx(i_3)==3
    E=2;
    F=2;
end
    s_5(i_3)=E*exp(j*2*pi*F/M);
    for ii_3=1:fs/srr
        s_6((i_3-1)*fs/srr+ii_3)=s_5(i_3)*exp(j*2*pi*((i_3-1)*fs/srr+ii_3)*fc2/fs);
    end
end
% ********************* signal4********************************
G=0;
H=0;
ffff=300;
N=ndd*fs/srr;
xxx=randint(nd,1,4);
M=2;
for i_4=1:nd
    if xxx(i_4)==0
    G=2;
    H=1/2;                               %如果将下面H的四个值改为0.5;1;1.5;2,即信号为QPSK 调制，这时在循环频率轴处峰值将消失；
elseif xxx(i_4)==1
    G=2;
    H=1/2;
elseif xxx(i_4)==2
    G=2;
    H=-1/2;
elseif xxx(i_4)==3
     G=2;
     H=-1/2;
    end
    s_7(i_4)=G*exp(j*2*pi*H*0);
    for ii_4=1:fs/sr
        s_8((i_4-1)*fs/sr+ii_4)=s_7(i_4)*exp(j*2*pi*((i_4-1)*fs/sr+ii_4)*(fc3+H*ffff)/fs);
    end
end
for i=1:length(s_2)
    s(i)=(s_2(i)+s_6(i));                       %合成OFDM信号
end
for i=1:length(s_2)
    s(i)=s(i);
end
sss=1/length(s)*abs(pwelch(s)).*abs(pwelch(s));
ss=ss+sss;
ss=1/100*ss;
% ********************* start ********************************
ff=25;                                           %搜索频率步进25Hz
     for tt=1:1:80                                %循环搜索频率点数
        r_1=0;
        r_2=0;
        for i=1:4000
            r_1=(r_1+s(i)*s(i)*conj(s(i))*exp(-j*2*pi*(1000+ff*tt)*i/fs));
            r_2=(r_2+s(i)*(s(i))*exp(-j*4*pi*(1000+ff*tt)*i/fs));
        end
        rr_1=1/4000*r_1;
        rr_2=1/4000*r_2;
        rr_3(tt)=rr_1;
        rr_4(tt)=rr_2;
     end
    rrr_1=abs((rr_3));
    shuju=rrr_1.*rrr_1;
    menxian=1/8*(shuju(28)+shuju(36)+shuju(44)+shuju(52));
    for iij=1:length(shuju)
        if shuju(iij)>menxian
            jieguo=jieguo+1;
        end
    end
    if jieguo>4
        jielun=jielun+1;
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    figure(1)
    subplot(2,1,1)
    f=1:2.5:length(ss)/2*2.5
    plot(f,ss(1:length(ss)/2).*ss(1:length(ss)/2)/10e4)
    xlabel('频率Hz')
    ylabel('幅值')
    title('重叠双信号频谱图')
    subplot(2,1,2)

    x=1:1:80
    plot(1000+x*25,(rrr_1.*rrr_1)/500000)      %平方的目的是为了幅值更明显
    xlabel('频率Hz')
    ylabel('幅值')
    title('不同频率处的三阶循环累积量值')  
%     colormap(1,1,1)
    figure(2)
    rrr_2=abs((rr_4))
    x=1:1:80
    plot(1000+x*25,rrr_2)      %平方的目的是为了幅值更明显
    xlabel('Hz')
    ylabel('幅值')
    title('循环频率轴截面图')  
    %%%%%%%%%%%%%%%%%%%%%%%%%%%
    q_1=0;
    q_2=0;
    q_3=0;
    q_4=0;
    q_5=0;
    q_6=0;    
%  a=1*fs/40:1*fs/40:12*fs/40
a_1=2*(1700);
a_2=2*(1900);
a_3=2*(1975);
a_4=2*(2025);
a_5=2*(2100);
a_6=2*(2300);
    for tt=1:800
        for i=1:4000
            q_1=q_1+s(i)*(s(i+tt))*exp(-j*2*pi*a_1*i/fs);
            q_2=q_2+s(i)*(s(i+tt))*exp(-j*2*pi*a_2*i/fs);
            q_3=q_3+s(i)*(s(i+tt))*exp(-j*2*pi*a_3*i/fs);
            q_4=q_4+s(i)*(s(i+tt))*exp(-j*2*pi*a_4*i/fs);
            q_5=q_5+s(i)*(s(i+tt))*exp(-j*2*pi*a_5*i/fs);
            q_6=q_6+s(i)*(s(i+tt))*exp(-j*2*pi*a_6*i/fs);       
        end
        qq_1(tt)=1/4000*q_1;
        qq_2(tt)=1/4000*q_2;
        qq_3(tt)=1/4000*q_3;
        qq_4(tt)=1/4000*q_4;
        qq_5(tt)=1/4000*q_5;
        qq_6(tt)=1/4000*q_6;
    end
    qqq_1=abs(fft(qq_1)).*abs(fft(qq_1))/100;
    qqq_2=abs(fft(qq_2)).*abs(fft(qq_2))/100;
    qqq_3=abs(fft(qq_3));
    qqq_4=abs(fft(qq_4));
    qqq_5=abs(fft(qq_5)).*abs(fft(qq_5))/100;
    qqq_6=abs(fft(qq_6)).*abs(fft(qq_6))/100;
    figure(1)
    subplot(3,2,1)
    plot((qqq_1(10:length(qqq_1))))
    xlabel('数据个数')
    ylabel('幅值')
    title('(a)频率1700Hz处的循环功率谱截面')  
    subplot(3,2,2)
    plot((qqq_2(10:length(qqq_2))))
    xlabel('数据个数')
    ylabel('幅值')
    title('(b)频率1900Hz处的循环功率谱截面')  
    subplot(3,2,3)
    plot((qqq_3(10:length(qqq_3))))
    xlabel('数据个数')
    ylabel('幅值')
    title('(c)频率1950Hz处的循环功率谱截面')  
    subplot(3,2,4)
    plot((qqq_4(10:length(qqq_4))))
    xlabel('数据个数')
    ylabel('幅值')
    title('(d)频率2050Hz处的循环功率谱截面')      
    subplot(3,2,5)
    plot((qqq_5(10:length(qqq_5))))
    xlabel('数据个数')
    ylabel('幅值')
    title('(e)频率2100Hz处的循环功率谱截面')      
    subplot(3,2,6)
    plot((qqq_6(10:length(qqq_6))))    
    xlabel('数据个数')
    ylabel('幅值')
    title('(f)频率2300Hz处的循环功率谱截面')