untitled_6.asv

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

clear
clc
sr=40       ;       % symbol rate没什么影响，不能设太大
fc=300;            %若与载波频率比值很大，则有较大影响
fs=2000    ;       % sample freq
Ts=1/fs;
nd=96       ;       % number of data没什么影响
srr=50        ;       % symbol rate
fc1=210    ;       % carrier freq1
fc2=310    ;       % carrier freq2     
% fc=(fc1+fc2)/2 ;
fss=2000    ;       % sample freq
ndd=120       ;       % number of data
loop=50        ;       % number of sample
taps=20        ;       % number of filter taps

% ********************* start ********************************

N=ndd*fss/srr;
xxx=randint(nd,1,4)
yy=dmod(xxx,fc,sr,fs,'fsk',4)
e_yy=mean(yy);
for i=1:ndd*fss/srr
yy(i)=(yy(i)-e_yy);
end
xx=randint(ndd,1,2 )
y=dmod(xx,fc2,srr,fss,'psk',2)
e_y=mean(y);
for i=1:ndd*fss/srr
y(i)=(y(i)-e_y);
end
for i=1:ndd*fss/srr
x(i)=(yy(i)+y(i));
end
% x=awgn(x,5)
r_1=0;
r_2=0;
r_3=0;
r_4=0;
r_5=0;
r_6=0;
%  a=1*fs/40:1*fs/40:12*fs/40
a_1=2*(fc2-1/2*sr);
a_2=2*(fc2);
a_3=2*(fc+1/2*sr);
a_4=2*(fc-1/2*sr);
a_5=2*(fc+3/2*sr);
a_6=2*(fc-3/2*sr);

    for tt=1:800
        for i=1:4000
            r_1=r_1+x(i)*(x(i+tt))*exp(j*2*pi*a_1*i/fs);
            r_2=r_2+x(i)*(x(i+tt))*exp(j*2*pi*a_2*i/fs);
            r_3=r_3+x(i)*(x(i+tt))*exp(j*2*pi*a_3*i/fs);
            r_4=r_4+x(i)*(x(i+tt))*exp(j*2*pi*a_4*i/fs);
            r_5=r_5+x(i)*(x(i+tt))*exp(j*2*pi*a_5*i/fs);
            r_6=r_6+x(i)*(x(i+tt))*exp(j*2*pi*a_6*i/fs);
        end
        rr_1(tt)=1/3200*r_1;
        rr_2(tt)=1/3200*r_2;
        rr_3(tt)=1/3200*r_3;
        rr_4(tt)=1/3200*r_4;
        rr_5(tt)=1/3200*r_5;
        rr_6(tt)=1/3200*r_6;        
    end
    rrr_1=abs(fft(rr_1)).*abs(fft(rr_1));
    rrr_2=abs(fft(rr_2)).*abs(fft(rr_2));
    rrr_3=abs(fft(rr_3)).*abs(fft(rr_3));
    rrr_4=abs(fft(rr_4)).*abs(fft(rr_4));
    rrr_5=abs(fft(rr_5)).*abs(fft(rr_5));
    rrr_6=abs(fft(rr_6)).*abs(fft(rr_6));    
% end
figure(1)
% tt=1:600
% a=1*fs/40:1*fs/40:12*fs/40
% mesh(a,tt,rrr)
subplot(3,2,1)
plot((rrr_1(10:length(rrr_1))))
title('在循环频率为载频的截面功率谱图') 
subplot(3,2,2)
plot((rrr_2(10:length(rrr_2))))
title('信号2在循环频率为载频的截面功率谱图') 
subplot(3,2,3)
plot((rrr_3(10:length(rrr_3))))
title('信号3在循环频率为载频的截面功率谱图') 
subplot(3,2,4)
plot((rrr_4(10:length(rrr_4))))
title('信号4在循环频率为载频的截面功率谱图') 
subplot(3,2,5)
plot((rrr_5(10:length(rrr_5))))
title('信号5在循环频率为载频的截面功率谱图') 
subplot(3,2,6)
plot((rrr_6(10:length(rrr_6))))
title('信号6在循环频率为载频的截面功率谱图') 
% figure(4)
% plot(abs(rr_2))
% hold on
% plot(abs(rr_2))