%變量說明，如下：
%x－待分析信號，以BPSK信號為例
%alpha－循環頻率
%fs－采樣頻率
%g－平滑窗函數，長度為平滑點數
%M－平滑窗長，即平滑點數
%L－一般設為L=M，即省略該參數的設置
%S-該切片的譜相關密度
%df-譜相關密度上的頻率分辨率



%初始化
clear;
clc;
%%%%%參量說明%%%%%%%%
pi=3.1415926;
fc=10;   %載波頻率
N=2000; %采樣點數
M=20;  %碼元數
L=N/M;  %每個碼元采樣點數
Rb=10;   %碼元速率
Tb=1/Rb;%碼元寬度
fs=40;
dt=1/fs;%采樣時域間隔
%采樣頻率
df=1/(dt*N);%采樣頻域間隔
t=[-(N*dt)/2+dt/2:dt:(N*dt)/2];
f=[-(N*df/2)+df/2:df:(N*df/2)];
%%%%%%NRZ碼%%%%%%%%%%%%%%%%%
a=floor(2*rand(1,M))*2-1;
for i=1:L
    s(i+[0:M-1]*L)=a;%采樣過程
end
%%%%%%%高斯預濾波%%%%%%%%%
Bb=Tb/0.3;%BT=0.3
alpha=sqrt(logm(2)/2/Bb^2);
H=exp(-alpha^2*f.^2);%高斯濾波器頻域表達式
S=fft(s);
S=[S(N/2+1:N),S(1:N/2)]*dt;
Y=S.*H;
Y=[Y(N/2+1:N),Y(1:N/2)];
y=ifft(Y)/dt;
y=real([y(N/2+1:N),y(1:N/2)]);
%%%%%%%串并轉換%%%%%%%%%%%
%It
for k=0:2*L:N-1;
kk=1:2:2*L;
kkk=1:L;
It(k+kk)=y(k+kkk+L);
It(k+kk+1)=y(k+kkk+L);
end
%Qt
for k=0:2*L:N-1;
kk=1:2:2*L;
kkk=1:L;
Qt(k+kk)=y(k+kkk);
Qt(k+kk+1)=y(k+kkk);
end
Itt=It.*cos(pi*t/2/Tb);
Qtt=Qt.*sin(pi*t/2/Tb);
gmsk=Itt.*cos(2*pi*fc*t)-Qtt.*sin(2*pi*fc*t);


%通過循環頻率的取值不同作出三維的循環譜圖
for alpha=0:0.1:3*fc

  
%設置平滑點數
N1=15; 
%設置平滑窗類型
g='hamming'; 
%設置L的值
if ~exist('L1') 
    		L1=N1;
end
if L1<1 
    error('L1 must be no less than unity.');
end
 
%為方便后續計算，將信號x從行向量變為列向量，并記錄相應行和列的個數
[row col]=size(gmsk);
 if col>2
    gmsk=gmsk';
    [row col]=size(gmsk);
 end


%所取時間段中的時間采樣的個數
lenx=row; 
%循環頻率的增量
d_alpha=fs/lenx; 
%頻率增量
d_f=fs/lenx;
%所設定的循環頻率對應的量化值
interval_f_N=round(alpha/d_alpha); %取最靠近的整數值
%計算平滑窗的個數
f_N=floor((lenx-interval_f_N-N1)/L1)+1;


%對信號x做fft變換
X=fftshift(fft(gmsk)); %將橫坐標零點放在中心(重要,幾乎每個程序都要這樣做)
%依次記錄平滑窗的編號，以便遍歷每個平滑窗
fnum=1:f_N; 
%為方便后續計算構造矩陣t，其列數為平滑點數，行數為平滑窗的個數
m=(1:N1)';
t1=zeros(N1,f_N);
t1=m(:,ones(f_N,1))+(fnum(ones(N1,1),:)-1)*L1;%從1開始每次取M個平滑點數

%對每個平滑窗采用非矩形窗函數，這個非矩形窗函數的長度為平滑點數
g=feval(g,N1); 
%為方便后續計算構造平滑窗函數矩陣
window_M=g(:,ones(f_N,1));
%x是一維情況，目前只用到這種情況
if col==1 
%初始化X1和X2
X1=zeros(N1,f_N);
X2=zeros(N1,f_N);
%計算X1
X1=X(t1).*window_M;
%計算X2
X2=X(t1+interval_f_N).*window_M; 
%x是多維情況
else 
    X1=X(t1,col);
    X1=reshape(X1,[N1,f_N]);
    X1=X1.*window_M;
    X2=X(t1+interval_f_N,1);
    X2=reshape(X2,[N1,f_N]);
    X2=X2.*window_M;
end
 
%不做平滑
if N1==1 
    S=conj(X1).*X2;
     S=S/lenx;
 
%做M點平滑
else
%計算譜相關
Stmp=conj(X1).*X2; 
%對平滑窗內求和
S=sum(Stmp); 
    S=S/lenx;     
end
 
%df是頻率增量
df=-1*(fs/2-d_alpha*floor(N1/2)-alpha/2):d_alpha*L1:(fs/2-d_alpha*floor(N1/2)-alpha/2);
%設定循環譜圖中df和S的范圍
df=df(1:min(length(S),length(df)));
S=S(1:min(length(S),length(df)));
%作圖
plot3(alpha*ones(1,length(df)),df,abs(S));
grid on;
hold on;
end
hold off;