Problems in doing STFT to GMSK signal

Hi all,

I'm trying to do STFT to GMSK signal with function "spectrogram.m", I can get the time-frequency plot but it seems not right. The problem is the center frequency seems not equal to the carrier frequency.

I am not sure whether i set the wrong parameter.

My code is listed below.

Thank you ALL!

Andrew

%CODE1:MAIN%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc
close all
clear

%%generate gmsk(gsm) signal

ini_phase=0*pi;
fc=900000000;%载波频率
sr=270833;  % Symbol rate调制速率
tb=1/sr;
nn=8;        %可调节采样率
ffs=sr*nn;       %sampling rate
nd = 1000;      % Number of symbols that simulates in each loop
ebn0=10;      % Eb/N0

s0 = gmsk_generate(nd,ffs,sr,fc,ebn0,ini_phase);
fs=ffs;
x=s0';
figure
[SS FF TT PP]=spectrogram(x,60,20,32,fs);   %how to set 'fs' ? Are the parameters suitable?

contour(TT*1e3,FF*1e-6,abs(SS));
axis tight;
view(0,90);
xlabel('times(ms)');ylabel('frequency(MHz)');title('时频图');
axis([0 3.5 0 2.2])

%%CODE2:generate gmask%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function s = gmsk_generate(nd,ffs,sr,fc,ebn0,ini_phase) %对符号数据升采样后再经过同样采样率的高斯滤波器得到MSK调制的输入数据 %GMSK为连续相位调制故再用该数据进行MSK的相位调制得到GMSK信号 %sr：符号速率 fc：载频 ffs：采样率

      mmfs=8;
    %-----根据采样率确定产生信号源时是否升采样
    %因为过采样后经过高斯滤波，若过采样过低不能较好的表示高斯滤波。
    %故若过采样过低则先用相对较高采样，最后再降采样
    if ffs/sr<8  
    fs=mmfs*ffs;
    else
        fs=ffs;
    end
    %-------------------------------------------
    ml=1;            % ml:Number of modulation levels 
    br=sr.*ml;       % Bit rate
    IPOINT=fs/sr;    % Number of oversamples

    %************************* Filter initialization ***************************

    irfn=21;                             % Number of taps
    Bb=0.3*sr;
    Bb2=0.6*sr;
    [xh] = gaussf(Bb,irfn,IPOINT,sr,1);   %Transmitter filter coefficients 
    [xh2] =gaussf(Bb2,irfn,IPOINT,sr,0);  %Receiver filter coefficients 

    %*************************** Data generation ******************************** 

        data=randint(nd,1,2);%产生基带信号
        x=2*data-1;
        nChan = size(x, 2);     % # of channels
        sigLen = size(x, 1);    % input signal length in symbols

     %------------------采样－－－－－－－－－－－－－－－－－－－－
             T0=0*rand/sr;      %采样初始时刻
            k1=fix(T0*fs);
            resT=T0-k1/fs;
             x1(1:k1)= x(1);

             if resT==0 && k1==0
            stt=1;
            else stt=2;        
             end
           ist=length(x1);
           restk=ist+1;
           restt=resT;
            for ii=stt:length(x)
                T=resT+1/sr;
                k1(ii)=fix(T*fs);
                resT=T-k1(ii)/fs;
                 x1(ist+1:ist+1+k1(ii)-1)= x(ii);
                 ist=ist+1+k1(ii)-1;
                 restk=[restk,ist+1];   %记录采样周期内包含两个符号的点的位置
                 restt=[restt,resT];    %记录采样周期内包含两个采样点与符号变化处的时间差
            end

    %*************************** GMSK Modulation ********************************  

        data1=x1;
        data2=conv(data1,xh)./IPOINT;                           % NEW for GMSK
        syncpoint =fix(irfn*IPOINT);

        th=zeros(1,length(data2)+1);
        ich2=zeros(1,length(data2)+1);
        qch2=zeros(1,length(data2)+1);

        k=1;
        for ii=2:length(data2)+1
            if ii==restk(k)+1+fix(syncpoint/2)    %判断是否在符号变化处的点
                if ii==2
                  th(1,ii)=th(1,ii-1)+pi/2*data2(1,ii-1)./IPOINT;
                  k=k+1;
                else
                  th(1,ii)=th(1,ii-1)+(1-restt(k)*fs)*pi/2*data2(1,ii-1)./IPOINT+restt(k)*fs*pi/2*data2(1,ii-2)./IPOINT;
                  k=k+1;
                  if k>length(restk)
                      k=length(restk);
                  end
                end
            else
            th(1,ii)=th(1,ii-1)+pi/2*data2(1,ii-1)./IPOINT;
            end
        end
        ich2=cos(th);
        qch2=sin(th);

    %************************** Attenuation Calculation ***********************

        spow=sum(ich2.*ich2+qch2.*qch2)/nd;  % sum: built in function
      attn=0.5*spow*sr/br*10.^(-ebn0/10);
      attn=sqrt(attn);                     % sqrt: built in function

    %********************* Add White Gaussian Noise (AWGN) **********************
      
        [ich3,qch3]= comb(ich2,qch2,attn);% add white gaussian noise
     % ich3=ich2;
      %qch3=qch2;

    %-------------相位检测------------    
        ich=ich2;
        qch=qch2;
        ich6=ich(fix(syncpoint/2)+1:length(ich)-fix(syncpoint/2));
        qch6=qch(fix(syncpoint/2)+1:length(qch)-fix(syncpoint/2));
        s = complex(ich6.', qch6.');
        gphase=angle(s);
        % figure(1)
      %plot(1:length(gphase),gphase)     %未降采样时的相位变化
      %grid on
      %hold on
      %axis([1,length(gphase),-pi,pi]) 

    if ffs/sr<8  
    mmfs=mmfs;
    else
        mmfs=1;
    end
      ich7=ich6(1:mmfs:end);
        qch7=qch6(1:mmfs:end);
        s2 = complex(ich7.', qch7.');
        gphase2=angle(s2);
         %figure(2)
     % plot(1:length(gphase2),gphase2,'r')    %降采样时后的相位变化
     % grid on
     % hold on
      %axis([1,length(gphase2),-pi,pi]) 
      %---------------------------------
       %-------接收滤波    
        [ich4,qch4] = compconv(ich3,qch3,xh2);
       %------------------
        syncpoint =fix(irfn*IPOINT);
        ich5=ich4(syncpoint:length(ich4)-syncpoint);
        qch5=qch4(syncpoint:length(qch4)-syncpoint);
       ss = complex(ich5.', qch5.');
       %-----根据采样率进行降采样
       if ffs/sr<8  
       s=ss(1:mmfs:end);
       else
        s=ss;
       end
       %------------------------

        kk=0:length(s)-1;
        s = s .* exp(j*ini_phase);    %加初始相位
        s=s.*exp(j*2*pi*fc*kk.'/fs);  %加载频  形成gmsk信号

      t=nd/sr;                       %仿真时间
      dt=t/length(s);                %仿真时间间隔
      s=s/max(s);                    %归一化