[Discuss-gnuradio] A GSM rough timing/frequency correction block.

[Achilleas Anastasopoulos]

Hi,

some time ago I got interested in receiving GSM with GnuRadio.

I attach a block that does rough timing  estimation and frequency 
correction based on FCCH blocks.
Just put it in your gr-howto directory and modify Makefile.am and 
howto.i accordingly.
Also attached is a python code that I used to test it.

I have not tested it with real data from the usrp yet (no DBSRX db in 
the lab), but with some online files (Robert's samples):
http://www.segfault.net/gsm/resources/

I would be interested to hear how the block performs with real-time
data from the usrp.

Next step (if/when I get some time) is to do fine timing estimation
jointly with channel estimation based on the Normal burst training sequence.

Best,
Achilleas

/* -*- c++ -*- */
/*
 * Copyright 2004 Free Software Foundation, Inc.
 *
 * This file is part of GNU Radio
 *
 * GNU Radio is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2, or (at your option)
 * any later version.
 *
 * GNU Radio is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with GNU Radio; see the file COPYING.  If not, write to
 * the Free Software Foundation, Inc., 51 Franklin Street,
 * Boston, MA 02110-1301, USA.
 */

#ifndef INCLUDED_CORRELATOR_STATE_TYPE_H
#define INCLUDED_CORRELATOR_STATE_TYPE_H

typedef enum {
  CORRELATOR_SEARCHING = 200, CORRELATOR_SLEEPING, CORRELATOR_TRACKING, 
CORRELATOR_SKIPPING
} correlator_state_type_t;

#endif

/* -*- c++ -*- */
/*
 * Copyright 2004 Free Software Foundation, Inc.
 * 
 * This file is part of GNU Radio
 * 
 * GNU Radio is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2, or (at your option)
 * any later version.
 * 
 * GNU Radio is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with GNU Radio; see the file COPYING.  If not, write to
 * the Free Software Foundation, Inc., 51 Franklin Street,
 * Boston, MA 02110-1301, USA.
 */

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include <howto_gsm_correlator_c.h>
#include <gr_io_signature.h>
#include <stdexcept>
#include <assert.h>
#include <iostream>
  

howto_gsm_correlator_c_sptr 
howto_make_gsm_correlator_c (
    int samples_per_symbol,
    int sequence_length,
    float threshold,
    int maxndata_past,
    int maxndata_future,
    int sleep_time,
    int max_track_time
    )
{
  return howto_gsm_correlator_c_sptr (
     new howto_gsm_correlator_c 
(samples_per_symbol,sequence_length,threshold,maxndata_past,maxndata_future,sleep_time,max_track_time)
         );
}

howto_gsm_correlator_c::howto_gsm_correlator_c (
    int samples_per_symbol,
    int sequence_length,
    float threshold,
    int maxndata_past,
    int maxndata_future,
    int sleep_time,
    int max_track_time
    )
  : gr_block ("gsm_correlator_c",
                          gr_make_io_signature (1, 1, sizeof (gr_complex)),
                          gr_make_io_signature (1, 1, sizeof (gr_complex))),  
  d_samples_per_symbol (samples_per_symbol),
  d_sequence_length (sequence_length),
  d_threshold(threshold),
  d_maxndata_past(maxndata_past),
  d_maxndata_future(maxndata_future),
  d_t(0),
  d_sleep_time(sleep_time),
  d_max_track_time(max_track_time),
  d_exponent(1.0),
  d_carrier(1.0)
{
    //set_output_multiple (1);
    //set_relative_rate (1.0); 
    set_history ( 
std::max(d_sequence_length*d_samples_per_symbol+1,(d_maxndata_past+d_max_track_time)*d_samples_per_symbol+1)
 );
    d_ndata = (d_maxndata_past+d_maxndata_future)*d_samples_per_symbol; // big 
number
    if (history()==1)
      d_state = CORRELATOR_SEARCHING;
    else
      d_state = CORRELATOR_SKIPPING;
    //printf("History is set to = %d\n\n",history());
}


void
howto_gsm_correlator_c::forecast (int noutput_items, gr_vector_int 
&ninput_items_required)
{
  int input_required =  noutput_items + history() - 1 ;
  ninput_items_required[0] = input_required;
  //printf("forecast %d %d\n",noutput_items,ninput_items_required[0]);
}


void gsm_correlate(const gr_complex *in, int Ls, float &mean, float &var) {
  mean = 0;
  for (int n=0;n<Ls;n++) {
    mean += arg(conj(in[-n-1])*in[-n]);
  }
  mean /= Ls;
  var = 0;
  float x;
  for (int n=0;n<Ls;n++) {
    x= fabs(arg(conj(in[-n-1])*in[-n])-mean);
    var += x*x; 
  }
  var/=Ls;
}



int
howto_gsm_correlator_c::general_work (int noutput_items,
                        gr_vector_int &ninput_items,
                        gr_vector_const_void_star &input_items,
                        gr_vector_void_star &output_items)
{
  float mean,var,omega;
  int L=history();
  //printf("history = %d\n",L);

  const gr_complex *in = (const gr_complex *) input_items[0];
  gr_complex *out = (gr_complex *) output_items[0];

  int nn=0;
  for (int n=0;n<noutput_items;n++) {

    // generate output
    if (d_ndata<(d_maxndata_past+d_maxndata_future)*d_samples_per_symbol) {
      assert(n+L-1-d_maxndata_past*d_samples_per_symbol-d_lag>=0);
      out[nn]=d_carrier*in[n+L-1-d_maxndata_past*d_samples_per_symbol-d_lag];
      nn++;
    }
    else {
      //out[nn]=0;
      //nn++;
      ;
    }

    d_time++;
    d_track_time++;
    d_ndata++; 
    d_carrier *= d_exponent;

    // correlate, track or sleep
    switch (d_state) {
      case CORRELATOR_SKIPPING:
        //printf("entered state CORRELATOR_SKIPPING\n");
        if(d_t>=L-2)
          d_state = CORRELATOR_SEARCHING;
      break;


      case CORRELATOR_SEARCHING:
        //printf("entered state CORRELATOR_SEARCHING\n");
        gsm_correlate(in+L-1+n,d_sequence_length*d_samples_per_symbol,mean,var);
        //printf("thresh=%f   mean=%f    var=%f\n",d_threshold,mean,var);
        if (var <= d_threshold) {
          if(d_max_track_time*d_samples_per_symbol>0) {
           d_state = CORRELATOR_TRACKING;
           d_min_var = var;
           d_time = 0;
           d_track_time = 0;
           omega = mean-M_PI/2/d_samples_per_symbol;
           d_exponent=gr_complex(cos(-omega),sin(-omega));
           printf("Found FCCH: %ld   %f %f \n",d_t,var,omega);
          }
          else {
           d_state = CORRELATOR_SLEEPING;
           d_rem_sleep_time=d_sleep_time*d_samples_per_symbol;
           d_ndata = 0;
           d_time = 0;
           d_lag = 0;
           //printf("LAG = 
%d,MAX_LAG=%d\n",d_lag,d_max_track_time*d_samples_per_symbol);
           assert(d_lag<=d_max_track_time*d_samples_per_symbol);
          }
        }
        break;


      case CORRELATOR_TRACKING:
        //printf("entered state CORRELATOR_TRACKING\n");
        gsm_correlate(in+L-1+n,d_sequence_length*d_samples_per_symbol,mean,var);
        //printf("thresh=%f   mean=%f var=%f  
track_time=%d\n",d_threshold,mean,var,d_track_time);
        if (var <= d_threshold && 
d_track_time<d_max_track_time*d_samples_per_symbol) {
           if (var < d_min_var) {
             d_min_var = var;
             d_time = 0;
             omega = mean-M_PI/2/d_samples_per_symbol;
             d_exponent=gr_complex(cos(-omega),sin(-omega));
             printf("Min variance update %ld    %d %f %f 
%f\n",d_t,d_track_time,var,d_min_var,omega);
           }
        }
        else {
           d_state = CORRELATOR_SLEEPING;  
           d_rem_sleep_time=d_sleep_time*d_samples_per_symbol-d_time;
           d_ndata = 0;
           d_lag = d_time; // at most d_max_track_time*d_samples_per_symbol-1
           //printf("LAG = 
%d,MAX_LAG=%d\n",d_lag,d_max_track_time*d_samples_per_symbol);
           assert(d_lag<=d_max_track_time*d_samples_per_symbol);
        }
        break;


      case CORRELATOR_SLEEPING:
        //printf("entered state CORRELATOR_SLEEPING\n");
        if (d_time >= d_rem_sleep_time)
           d_state = CORRELATOR_SEARCHING;
        break;

      default:
        throw std::runtime_error ("Invalid state type.");

    } // end switch

    d_t++;

  }
  

  consume_each(noutput_items);

  //printf("%d  %d\n",nn,noutput_items);
  return nn;
  //return noutput_items;
}

/* -*- c++ -*- */
/*
 * Copyright 2004 Free Software Foundation, Inc.
 * 
 * This file is part of GNU Radio
 * 
 * GNU Radio is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2, or (at your option)
 * any later version.
 * 
 * GNU Radio is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with GNU Radio; see the file COPYING.  If not, write to
 * the Free Software Foundation, Inc., 51 Franklin Street,
 * Boston, MA 02110-1301, USA.
 */

#ifndef INCLUDED_HOWTO_GSM_CORRELATOR_C_H
#define INCLUDED_HOWTO_GSM_CORRELATOR_C_H

#include <vector>
#include <gr_complex.h>
#include <gr_block.h>
#include <correlator_state_type.h>

class howto_gsm_correlator_c;
typedef boost::shared_ptr<howto_gsm_correlator_c> howto_gsm_correlator_c_sptr;

howto_gsm_correlator_c_sptr howto_make_gsm_correlator_c (
    int samples_per_symbol,
    int sequence_length,
    float threshold,
    int maxndata_past,
    int maxndata_future,
    int sleep_time,
    int max_track_time
);



class howto_gsm_correlator_c : public gr_block
{
  int d_samples_per_symbol;
  int d_sequence_length;
  float d_threshold;
  int d_maxndata_past;
  int d_maxndata_future;
  int d_sleep_time;
  int d_max_track_time;

  correlator_state_type_t d_state;
  int d_time;
  int d_rem_sleep_time;
  int d_track_time;
  int d_lag;
  long int d_t;
  int d_ndata;
  float d_min_var;
  gr_complex d_exponent;
  gr_complex d_carrier;

  friend howto_gsm_correlator_c_sptr howto_make_gsm_correlator_c (
    int samples_per_symbol,
    int sequence_length,
    float threshold,
    int maxndata_past,
    int maxndata_future,
    int sleep_time,
    int max_track_time
  );


  howto_gsm_correlator_c (
    int samples_per_symbol,
    int sequence_length,
    float threshold,
    int maxndata_past,
    int maxndata_future,
    int sleep_time,
    int max_track_time
  );


public:
  int samples_per_symbol() const { return d_samples_per_symbol; }
  int sequence_length () const { return d_sequence_length; }
  float threshold() const { return d_threshold; }
  int maxndata_past() const { return d_maxndata_past; }
  int maxndata_future() const { return d_maxndata_future; }
  int sleep_time() const { return d_sleep_time; }
  int max_track_time() const { return d_max_track_time; }

  void forecast (int noutput_items, gr_vector_int &ninput_items_required);

  int general_work (int noutput_items,
                    gr_vector_int &ninput_items, 
                    gr_vector_const_void_star &input_items,
                    gr_vector_void_star &output_items);
};

#endif

/* -*- c++ -*- */
/*
 * Copyright 2004 Free Software Foundation, Inc.
 *
 * This file is part of GNU Radio
 *
 * GNU Radio is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2, or (at your option)
 * any later version.
 *
 * GNU Radio is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with GNU Radio; see the file COPYING.  If not, write to
 * the Free Software Foundation, Inc., 51 Franklin Street,
 * Boston, MA 02110-1301, USA.
 */

GR_SWIG_BLOCK_MAGIC(howto,gsm_correlator_c);

howto_gsm_correlator_c_sptr howto_make_gsm_correlator_c (
    int samples_per_symbol,
    int sequence_length,
    float threshold,
    int maxndata_past,
    int maxndata_future,
    int sleep_time,
    int max_track_time
);


class howto_gsm_correlator_c : public gr_block
{
private:
  gsm_correlator_c (
    int samples_per_symbol,
    int sequence_length,
    float threshold,
    int maxndata_past,
    int maxndata_future,
    int sleep_time,
    int max_track_time
  );

public:
    int samples_per_symbol() const { return d_samples_per_symbol; }
    int sequence_length () const { return d_sequence_length; }
    float threshold () const { return d_threshold; }
    int maxndata_past() const { return d_maxndata_past; }
    int maxndata_future() const { return d_maxndata_future; }
    int sleep_time() const { return d_sleep_time; }
    int max_track_time() const { return d_max_track_time; }
};

#!/usr/bin/env python
from gnuradio import gr, gru, blks, optfir
from Numeric import *
from gnuradio.wxgui import slider, powermate
from gnuradio.wxgui import stdgui, fftsink, form, scopesink
import wx
from gnuradio import howto

class my_graph(gr.flow_graph):
    def __init__(self):
        gr.flow_graph.__init__(self)

#class my_graph (stdgui.gui_flow_graph):
    #def __init__(self,frame,panel,vbox,argv):
        #stdgui.gui_flow_graph.__init__ (self,frame,panel,vbox,argv)

        self.src = 
gr.file_source(gr.sizeof_gr_complex,'GSMSP_20070204_robert_dbsrx_953.6MHz_64.cfile',False);
        decimation = 64
        fs_a2d = 64e6
        fs_in=fs_a2d/decimation
        BW=200e3
        OSR = 4 # make sampling rate = 4 x Rb. 
        #Recall: Tb = (15/26) msec / 156.25 bits = 3.6923e-006 (approx)
        # so Ts = Tb/4 = 9.2308e-007 (approx) = 12/13 1e-6 (approx)
        a=13*OSR*decimation
        b=12*4*64
        c=gru.gcd(a,b)
        fs=fs_in*a/b

        chan_filt_coeffs = optfir.low_pass (1,           # gain
                                            fs_in,   # sampling rate
                                            0.8*BW,        # passband cutoff
                                            1.0*BW,       # stopband cutoff
                                            0.1,         # passband ripple
                                            60)          # stopband attenuation
        print "Number of coeffs = ", len(chan_filt_coeffs)
        self.filter = gr.fir_filter_ccf (1, chan_filt_coeffs)
        self.resample = blks.rational_resampler_ccf(self, a/c, b/c)
        self.agc = gr.agc_cc(1e-6,1.0,1.0,0.0);

        if 0:
          self.fft1 = fftsink.fft_sink_c (self, panel, 
                title = 'FFT 1', 
                fft_size=512, 
                y_per_div=10, ref_level=60, 
                fft_rate=5, 
                average=True, avg_alpha=0.01, 
                sample_rate=fs_in)
          vbox.Add (self.fft1.win, 10, wx.EXPAND)

          self.fft2 = fftsink.fft_sink_c (self, panel, 
                title = 'FFT 2', 
                fft_size=512, 
                y_per_div=10, ref_level=0, 
                fft_rate=5, 
                average=True, avg_alpha=0.01, 
                sample_rate=fs)
          vbox.Add (self.fft2.win, 10, wx.EXPAND)


        self.corr = howto.gsm_correlator_c(OSR, # samples per symbol
                142, # number of symbols for freq estimation
                0.01, # threshold
                8+3+142,# past symbols to output
                int(round(3+0.25+156.25*79)), # future symbols to output
                int(round(3+8.25+156.25*79))-200-50, # sleep time (in symbols)
                50) # track time (in symbols)
        #self.c2mag =gr.complex_to_mag()
        #self.scope = scopesink.scope_sink_f(self, panel, 
                #sample_rate=Q,
                #frame_decim=10,
                #v_scale=100,
                #t_scale=options.t_scale
         #)
        #vbox.Add (self.scope.win, 10, wx.EXPAND)
        self.sink = gr.file_sink(gr.sizeof_gr_complex,'gsm_64.cfile')

        

        # Connect  everything together
        self.connect 
(self.src,self.filter,self.resample,self.agc,self.corr,self.sink)
        #self.connect (self.src,self.fft1)
        #self.connect (self.agc,self.fft2)



if __name__ == '__main__':
    #app = stdgui.stdapp (my_graph, "Transmitter")
    #app.MainLoop ()  
    fg=my_graph()
    fg.run()