gsm_cr.m

自己编写的GMSK信号的循环谱密度程序

%%%%%%%%%%%%%%%%%%%%%GMSK%%%%%%%%%%%%%%%%%%%%
clear all;
%%%%%参量说明%%%%%%%%
pi=3.1415926;
fc=10;   %载波频率  10MHZ
N=2000; %采样点数
M=500;  %码元数
L=N/M;  %每个码元采样点数
Rb=4;   %码元速率  4MHZ/S
Tb=1/Rb;%码元宽度
dt=Tb/L;%采样时域间隔
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);
gmsk=normrnd(0,1,1,length(gmsk));

Npnt=128;
N1=100;
Xtnf=zeros((N-N1+1),Npnt);
for t=1:1:(N-N1+1)
    %function datafft= WinFFT(data,Npntfft)
    
    Xtnf(t,:)=WinFFT(gmsk(t:(t-1+N1)),Npnt);
end

Nmax=floor((Npnt-1)/2);
alphaN=-2*Nmax:2:2*Nmax;
fN=-Nmax:1:Nmax;
XX=zeros(length(fN),length(alphaN));

for t=1:1:(N-N1+1)
    for alpha=alphaN        
        for f=-(Nmax-abs(alpha/2)):1:(Nmax-abs(alpha/2))             
            if f+alpha/2 >= 0               
                x1=Xtnf(t,f+alpha/2+1);               
            else
                x1=Xtnf(t,Npnt+(f+alpha/2+1));
            end
            if f-alpha/2 >= 0
                x2=Xtnf(t,f-alpha/2+1);
            else               
                x2=Xtnf(t,Npnt+(f-alpha/2+1));
            end
            XX(f+Nmax+1,alpha/2+Nmax+1)=XX(f+Nmax+1,alpha/2+Nmax+1)+abs(x1*conj(x2));
            %if f==0
            %    if abs(x1)~=abs(x2)
            %        1111
            %    end
            %end          
        end
    end
   process_percentage=t/(N-N1+1)*100
end
XX=XX/((N-N1+1));
surfc(fN,alphaN,abs(XX(fN+Nmax+1,alphaN/2+Nmax+1)));