gsm_cr_2.m

待分析信号

%变量说明，如下：
%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;