常规信号的循环谱(代码部分)

前不久团队贴出一些常用信号的循环谱仿真图（见【GreenSim原创】常规信号的循环谱），由于篇幅所限这个帖子把剩余的几幅图和程序贴出来，供大家交流。

64QAM信号循环谱图

16QAM信号循环谱图

256QAM信号循环谱图
计算信号循环谱的函数：
这个函数用来计算数字调制信号的循环谱图，结果存放在sxf中，surf和shading interp这两个函数配合使用就可以得到信号的循环谱图，M一般取64，N取2048。需要注意的是，数字信号要保留一定的载频，不能是基带信号。
function sxf=MyCS_density(s,M,N)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  谱相关的计算函数
%%  参数说明：M是平滑窗,必须为2的幂
%%            N是信号序列长度,必须为2的幂
%%            s是时域序列
%%            sxf是返回参数
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
spe_s = fft(s,N);
k = N/2-M;
sxf = zeros(2*k+1,4*k+1);                          %因为菱形区之外的数为零
for ff = -k:k;
for alpha= -k+ff:1:2* ff                        %（菱形区）
for v = -M:M
temp1 = ff+v+alpha/2;
temp2 = ff+v-alpha/2;
if(temp1<1)
temp1 = temp1+N;
end
if(temp2<1)
temp2 =temp2+N;
end
cospe(v+M+1) = spe_s(temp1)* spe_s(temp2);
end
j=ff+k+1;
u=alpha+k+1;
sxf(j,u) =(abs( cospe))/(2*M*N);
end
end
计算信号循环谱alpha截面的函数：
这个函数用来计算数字调制信号的alpha截面图，结果存放在sxf_amp中，参数的设置和上一个函数一样，需要注意的是ff参数必须是数字频率。这个函数非常简洁也非常有用，尤其是在设计基于循环谱的参数估计和调制识别算法中，一般用这个函数。
function sxf_amp = MyCSalpha_imag_density(s,M,N,ff)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  谱相关alpha截面图的计算函数
%%  参数说明：M是平滑窗
%%            N是信号序列长度
%%            s是信号序列
%%            sxf_amp是幅度谱
%%            ff是频率
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
spe_s = fft(s,N);
spe_s=spe_s./2;
k = N/2-M;
sxf_amp = zeros(1,2*k+1);
for alpha=0:1:k-abs(ff))
for v =-M:M
temp1 = ff+v+alpha/2;
temp2 = ff+v-alpha/2;
if(temp1<1)
temp1 = abs(temp1);
A(temp1+1)=conj(spe_s(temp1+1));
else
A(temp1+1)=spe_s(temp1+1);
end
if(temp2<1)
temp2 = abs(temp2);
B(temp2+1)=conj(spe_s(temp2+1));
else
B(temp2+1)=spe_s(temp2+1);
end
cospe(v+M+1) = A(temp1+1)*B(temp2+1);
end
u=alpha+1;
sxf_amp(u) =(abs(cospe))/(2*M*N);
end
计算信号循环谱f截面的函数：
这个函数和上一个函数用法差不多，主要在设计参数估计和调制识别算法中使用。
function sxf_amp = MyCSf_imag_density(s,M,N,alpha)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  谱相关f截面图的计算函数
%%  参数说明：M是平滑窗
%%            N是信号序列长度
%%            s是复数序列
%%            sxf_amp是幅度谱
%%            alpha是数字频率
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
spe_s = fft(s,N);
spe_s=spe_s./2;
k = N/2-M;
sxf_amp = zeros(k+1,1);
for ff = 0:k-abs(alpha/2);
for v =-M:M
temp1 = ff+v+ alpha/2;
temp2 = ff+v- alpha/2;
if(temp1<1)
temp1 = abs(temp1);
A(temp1+1)=conj(spe_s(temp1+1));
else
A(temp1+1)=spe_s(temp1+1);
end
if(temp2<1)
temp2 = abs(temp2);
B(temp2+1)=conj(spe_s(temp2+1));
else
B(temp2+1)=spe_s(temp2+1);
end
cospe(v+M+1) = A(temp1+1)* B(temp2+1);
end
j=ff+1;
sxf_amp(j) =(abs (cospe))/(2*M*N);
end