clear;
clc;
fc=30*10^6;%載頻
rb=10*10^6;%碼元速率
t0=1*10^(-6);
fs=160*10^6;
n=fs/rb;%%%%%%一個碼元的采樣數
num=0;
%for k=1:10
xb=randint(1,100);
x2=xb.*2-1;%%%%%%%%%%%（0，1）—>(-1,1) 
t=1:length(x2)*n;
%%%%%%%%%%數據采樣
for i=1:length(x2)
    for j=1:n
        x3((i-1)*n+j)=x2(i);
    end
end
sig=cos(2*pi*fc*t/fs+x3.*pi/2);
x=sig(1:1000);
%x=awgn(x,10,'measured');
%%%%%%%%%幅度譜
%x=hilbert(x1);
[S,rho]=cyclic_cross_periodogram(x,x,64,'hamming','hamming',1);
N=64;
f=linspace(-0.5,0.5,N+1);
alpha=linspace(-0.5,0.5,N+1);
mesh(alpha,f,S)
xlabel('alpha/Hz');
ylabel('f/Hz');
zlabel('Amplitude');
S1=rho(floor(N/2),:)';
m1=max(rho*S1);
S2=S1/m1;%信號的歸一化循環譜密度
[S3,index1]=max(S2);%提取S2的最大峰值
m=mean(S2);%計算S2的均值
for k=1:N
    S4(k)=k*S2(k);
end
S5=sum(S4);
S6=sum(S2);
index2=S5/S6;%S2的中心位置
d=abs(index2-index1);%index1是S2峰值位置，d為S2的中心位置和峰值位置的差的絕對值
p=S2;%alpha軸上的周期譜線
s=sum(p.*p)/N;%Sx在alpha軸上的周期譜線平均能量S
v=sqrt(s-(sum(p)/N)^2);%標準偏差
k=mean(S2.^3)/v^3;%S2的偏度系數
for i=1:N
    for j=1:N
if rho(i,j)==1
    rho(i,j)=0;
end
    end
end
[rho1,index3]=max(max(rho));%rho的次大值
if rho1<0.99&r3>0.15&Smax>0.5&m>0.1
    num=num+1;
end
%end