fskc31.asv

在工程里识别不同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
for xuanhuan=1:100
% ********************* 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)+2*s_6(i));                       %合成OFDM信号
end
for i=1:length(s_2)
    s(i)=s(i);
end
sss=abs(fft(s));
ss=sss+ss;
end
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
    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))
    xlabel('频率Hz')
    ylabel('幅值')
    title('重叠双信号频谱图')
    subplot(2,1,2)
    rrr_1=abs((rr_3))
    x=1:1:80
    plot(1000+x*25,(rrr_1.*rrr_1)/500)      %平方的目的是为了幅值更明显
    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('循环频率轴截面图')