Added new blocks written in python for new experimental gsm receiver.

FCCH burst tagger is element of hierarchical block - FCCH detector that adds tags after every detected FCCH burst. The value of each tag is a frequency offset estimate.
SCH detector waits for tags from FCCH detector which are used to find SCH burst position. It is unfinished.

5 files changed, 336 insertions, 1 deletions

diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt
index 8640d73..2c76806 100644
--- a/python/CMakeLists.txt
+++ b/python/CMakeLists.txt
@@ -31,7 +31,10 @@ endif()
 GR_PYTHON_INSTALL(
     FILES
     __init__.py
-    receiver_hier.py DESTINATION ${GR_PYTHON_DIR}/gsm
+    receiver_hier.py
+    fcch_burst_tagger.py
+    sch_detector.py
+    fcch_detector.py DESTINATION ${GR_PYTHON_DIR}/gsm
 )
 
 ########################################################################
diff --git a/python/__init__.py b/python/__init__.py
index d6dec38..ee72a0b 100644
--- a/python/__init__.py
+++ b/python/__init__.py
@@ -46,6 +46,10 @@ from gsm_swig import *
 
 # import any pure python here
 from receiver_hier import receiver_hier
+
+from fcch_burst_tagger import fcch_burst_tagger
+from sch_detector import sch_detector
+from fcch_detector import fcch_detector
 #
 
 # ----------------------------------------------------------------
diff --git a/python/fcch_burst_tagger.py b/python/fcch_burst_tagger.py
new file mode 100644
index 0000000..5142096
--- /dev/null
+++ b/python/fcch_burst_tagger.py
@@ -0,0 +1,125 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+# 
+# Copyright 2014 Piotr Krysik pkrysik@elka.pw.edu.pl
+# 
+# This 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 3, or (at your option)
+# any later version.
+# 
+# This software 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 this software; see the file COPYING.  If not, write to
+# the Free Software Foundation, Inc., 51 Franklin Street,
+# Boston, MA 02110-1301, USA.
+# 
+
+from numpy import *
+from pylab import *
+from gnuradio import gr
+import pmt
+from scipy.signal.chirpz import ZoomFFT
+
+class fcch_burst_tagger(gr.sync_block):
+    """
+    docstring for block fcch_burst_tagger
+    """
+    def __init__(self, OSR):
+        gr.sync_block.__init__(self,
+            name="fcch_burst_tagger",
+            in_sig=[complex64, float32],
+            out_sig=[complex64])
+
+        self.state=False
+        self.symbol_rate = 1625000/6
+        self.OSR=OSR
+        self.samp_rate = self.symbol_rate*OSR
+        self.burst_size = int(156.25*self.OSR)
+        self.guard_period = int(round(8.25*self.OSR))
+        self.block_size = self.burst_size+self.guard_period
+        self.processed_block_size = int(142*self.OSR)
+        self.set_history(self.burst_size)
+        self.set_output_multiple(self.guard_period)
+        self.prev_offset=0
+        
+        #parameters of zoomfft frequency estimator
+        f1 = self.symbol_rate/4*0.9
+        f2 = self.symbol_rate/4*1.1
+        m=5000*self.OSR
+        self.zoomfft = ZoomFFT(self.processed_block_size, f1, f2, m, Fs=self.samp_rate)
+        self.f_axis = linspace(f1,f2,m)
+
+    def work(self, input_items, output_items):
+        in0=input_items[0]
+        output_items[0][:] = in0[self.history()-1:]
+
+        threshold = input_items[1][self.history()-1:]
+        threshold_diff = diff(concatenate([[0],threshold]))
+        up_to_high_indexes = nonzero(threshold_diff>0)[0]
+
+        up_to_high_idx=[] 
+        
+        for up_to_high_idx in up_to_high_indexes:           #look for "high" value at the trigger
+            if up_to_high_idx==0 and self.state==True:    #if it's not transition from "low" to "high"
+                continue                                    #then continue
+            self.state=True                               #if found - change state
+        
+        if self.state==True and up_to_high_idx and any(threshold_diff<0):          #and look for transition from high to low
+            last_up_to_high_idx = up_to_high_idx
+            last_high_to_low_idx = nonzero(threshold_diff<0)[0][-1]
+            
+            if last_high_to_low_idx-last_up_to_high_idx>0:
+                coarse_idx = int(last_high_to_low_idx+self.history()-self.guard_period-self.burst_size)
+                inst_freq = angle(in0[coarse_idx:coarse_idx+self.block_size]*in0[coarse_idx-self.OSR:coarse_idx+self.block_size-self.OSR].conj())/(2*pi)*self.symbol_rate #instantaneus frequency estimate
+                precise_idx = self.find_best_position(inst_freq)
+#                measured_freq = mean(inst_freq[precise_idx:precise_idx+self.processed_block_size])
+                expected_freq = self.symbol_rate/4
+                zoomed_spectrum = abs(self.zoomfft(in0[coarse_idx+precise_idx:coarse_idx+precise_idx+self.processed_block_size]))
+                measured_freq = self.f_axis[argmax(zoomed_spectrum)]
+                freq_offset = measured_freq - expected_freq
+                offset = self.nitems_written(0) + coarse_idx + precise_idx - self.guard_period
+                key = pmt.string_to_symbol("fcch")
+                value =  pmt.from_double(freq_offset)
+                self.add_item_tag(0,offset, key, value)
+                self.state=False
+
+#   Some additional plots and prints for debugging
+#                print "coarse_idx+precise_idx",coarse_idx+precise_idx
+#                print "offset-self.nitems_written(0):",offset-self.nitems_written(0)
+                print offset-self.prev_offset
+                self.prev_offset=offset
+                print "freq offset", freq_offset
+#                freq_offset = measured_freq - expected_freq
+#                plot(self.f_axis, zoomed_spectrum)
+#                show()
+#                plot(inst_freq[precise_idx:precise_idx+self.burst_size])
+#                show()
+#                plot(unwrap(angle(in0[coarse_idx+precise_idx:coarse_idx+precise_idx+self.burst_size])))
+#                show()
+#                
+        return len(output_items[0])
+
+    def find_best_position(self, inst_freq):
+        lowest_max_min_diff = 1e6 #1e6 - just some large value
+        start_pos = 0
+        
+        for ii in xrange(0,int(30*self.OSR)):
+            min_inst_freq = min(inst_freq[ii:self.processed_block_size+ii-1]);
+            max_inst_freq = max(inst_freq[ii:self.processed_block_size+ii-1]);
+
+            if (lowest_max_min_diff > max_inst_freq - min_inst_freq):
+                lowest_max_min_diff = max_inst_freq - min_inst_freq;
+                start_pos = ii
+#                print 'start_pos',start_pos
+        
+#        plot(xrange(start_pos,start_pos+self.processed_block_size),inst_freq[start_pos:start_pos+self.processed_block_size],'r.')
+#        hold(True)
+#        plot(inst_freq)
+#        show()
+        
+        return start_pos
diff --git a/python/fcch_detector.py b/python/fcch_detector.py
new file mode 100644
index 0000000..627dd00
--- /dev/null
+++ b/python/fcch_detector.py
@@ -0,0 +1,71 @@
+#!/usr/bin/env python
+##################################################
+# Gnuradio Python Flow Graph
+# Title: FCCH Bursts Detector
+# Author: Piotr Krysik
+#
+# Description: Detects positions of FCCH bursts. At the end of each 
+# detected FCCH burst adds to the stream a tag with key "fcch" and value 
+# which is a frequency offset estimate. The input sampling frequency 
+# should be integer multiply of GSM GMKS symbol rate - 1625000/6 Hz.
+##################################################
+
+from gnuradio import blocks
+from gnuradio import gr
+from gnuradio.filter import firdes
+import gsm
+
+class fcch_detector(gr.hier_block2):
+
+    def __init__(self, OSR=4):
+        gr.hier_block2.__init__(
+            self, "FCCH bursts detector",
+            gr.io_signature(1, 1, gr.sizeof_gr_complex*1),
+            gr.io_signature(1, 1, gr.sizeof_gr_complex*1),
+        )
+
+        ##################################################
+        # Parameters
+        ##################################################
+        self.OSR = OSR
+
+        ##################################################
+        # Variables
+        ##################################################
+        self.f_symb = f_symb = 1625000.0/6.0
+        self.samp_rate = samp_rate = f_symb*OSR
+
+        ##################################################
+        # Blocks
+        ##################################################
+        self.gsm_fcch_burst_tagger_0 = gsm.fcch_burst_tagger(OSR)
+        self.blocks_threshold_ff_0_0 = blocks.threshold_ff(0, 0, 0)
+        self.blocks_threshold_ff_0 = blocks.threshold_ff(int((138)*samp_rate/f_symb), int((138)*samp_rate/f_symb), 0)
+        self.blocks_multiply_conjugate_cc_0 = blocks.multiply_conjugate_cc(1)
+        self.blocks_moving_average_xx_0 = blocks.moving_average_ff(int((142)*samp_rate/f_symb), 1, int(1e6))
+        self.blocks_delay_0 = blocks.delay(gr.sizeof_gr_complex*1, int(OSR))
+        self.blocks_complex_to_arg_0 = blocks.complex_to_arg(1)
+
+        ##################################################
+        # Connections
+        ##################################################
+        self.connect((self, 0), (self.blocks_multiply_conjugate_cc_0, 0))
+        self.connect((self.blocks_delay_0, 0), (self.blocks_multiply_conjugate_cc_0, 1))
+        self.connect((self.blocks_complex_to_arg_0, 0), (self.blocks_threshold_ff_0_0, 0))
+        self.connect((self, 0), (self.blocks_delay_0, 0))
+        self.connect((self.blocks_multiply_conjugate_cc_0, 0), (self.blocks_complex_to_arg_0, 0))
+        self.connect((self.blocks_moving_average_xx_0, 0), (self.blocks_threshold_ff_0, 0))
+        self.connect((self.blocks_threshold_ff_0_0, 0), (self.blocks_moving_average_xx_0, 0))
+        self.connect((self.gsm_fcch_burst_tagger_0, 0), (self, 0))
+        self.connect((self, 0), (self.gsm_fcch_burst_tagger_0, 0))
+        self.connect((self.blocks_threshold_ff_0, 0), (self.gsm_fcch_burst_tagger_0, 1))
+
+    def get_OSR(self):
+        return self.OSR
+
+    def set_OSR(self, OSR):
+        self.OSR = OSR
+        self.set_samp_rate(self.f_symb*self.OSR)
+        self.blocks_delay_0.set_dly(int(self.OSR))
+
+
diff --git a/python/sch_detector.py b/python/sch_detector.py
new file mode 100644
index 0000000..d80d032
--- /dev/null
+++ b/python/sch_detector.py
@@ -0,0 +1,132 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+# 
+# Copyright 2014 <+YOU OR YOUR COMPANY+>.
+# 
+# This 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 3, or (at your option)
+# any later version.
+# 
+# This software 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 this software; see the file COPYING.  If not, write to
+# the Free Software Foundation, Inc., 51 Franklin Street,
+# Boston, MA 02110-1301, USA.
+# 
+
+from numpy import *
+from pylab import *
+from gnuradio import gr
+import pmt
+from scipy.ndimage.filters import uniform_filter1d
+
+class sch_receiver():
+    """
+    docstring for class sch_reciever
+    """
+    def __init__(self, OSR):
+        self.sync_seq = array([1, 0, 1, 1, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0,
+                               0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1,
+                               0, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1,
+                               0, 1, 1, 1, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 1, 1])
+        self.OSR = OSR
+        sync_seq_msk_tmp = self.msk_mod(self.sync_seq, -1j)
+        self.sync_seq_msk = sync_seq_msk_tmp[5:59]
+        self.sync_seq_msk_interp = zeros(self.OSR*len(self.sync_seq_msk), dtype=np.complex64)
+        self.sync_seq_msk_interp[::OSR] = self.sync_seq_msk
+        self.L = 5
+
+    def msk_mod(self, x, start_point):
+        x_nrz = 2*x-1 
+        x_diffenc = x_nrz[1:]*x_nrz[0:-1]
+        mod_tmp = concatenate((array([start_point]),1j*x_diffenc))
+        return cumprod(mod_tmp)
+    
+    def get_chan_imp_resp(self, sch_burst):
+        sch_burst_bl = resize(array(sch_burst), (int(len(sch_burst)/self.OSR),self.OSR))
+        correlation_bl = zeros(shape(sch_burst_bl), dtype=np.complex64)
+        for ii in xrange(0,self.OSR):
+            correlation_bl[:,ii]=correlate(sch_burst_bl[:,ii],self.sync_seq_msk,'same')
+        
+        correlation_bl = correlation_bl/len(self.sync_seq_msk)
+        power_bl_mov_avg = uniform_filter1d(abs(correlation_bl)**2,5,mode='constant',axis=0)
+
+        print "correlation_bl.argmax()",argmax(abs(hstack(correlation_bl)))
+        print "power_bl_mov_avg.argmax()",hstack(power_bl_mov_avg).argmax()
+        print 'unravel_index(correlation_bl.argmax(), correlation_bl.shape)',unravel_index(correlation_bl.argmax(), correlation_bl.shape)
+        print 'unravel_index(power_bl_mov_avg.argmax(), power_bl_mov_avg.shape)',unravel_index(power_bl_mov_avg.argmax(), power_bl_mov_avg.shape)
+#        correlation = zeros(shape(sch_burst))
+#        correlation = correlate(sch_burst, self.sync_seq_msk_interp,'same')/len(self.sync_seq_msk)
+#        print "pozycja maksimum",argmax(abs(correlation))
+        plot(abs(hstack(correlation_bl))*1000)
+#        hold(True)
+#        plot(abs(sch_burst)*500)
+#        print shape(range(0,len(sch_burst),self.OSR))
+#        print shape(correlation_bl[:,0])
+        for ii in range(0,self.OSR):
+            plot(range(ii,len(correlation_bl[:,0])*self.OSR,self.OSR),power_bl_mov_avg[:,ii]*5e6,'r.')
+        show()
+    def receive(self, input_corr, chan_imp_resp):
+        pass
+
+class sch_detector(gr.sync_block):
+    """
+    docstring for block sch_detector
+    """
+    def __init__(self, OSR):
+        gr.sync_block.__init__(self,
+            name="sch_detector",
+            in_sig=[complex64],
+            out_sig=[complex64])
+        self.OSR = OSR
+        self.states = {"waiting_for_fcch_tag":1, "reaching_sch_burst":2, "sch_at_input_buffer":3}
+        self.state = self.states["waiting_for_fcch_tag"]
+        self.sch_offset = -100 #-100 - just some invalid value of sch burst position in the stream
+        self.burst_size = int(round(156.25*self.OSR))
+        self.guard_period = int(round(8.25*self.OSR))
+        self.block_size = self.burst_size + self.guard_period
+        self.set_history(self.block_size)
+        self.set_output_multiple(self.guard_period)
+        self.sch_receiver = sch_receiver(OSR)
+        
+    def work(self, input_items, output_items):
+        in0 = input_items[0]
+        out = output_items[0]
+        to_consume = len(in0)-self.history()
+        
+        if self.state == self.states["waiting_for_fcch_tag"]:
+            fcch_tags = []
+            
+            start = self.nitems_written(0)
+            stop = start + len(in0)
+            key = pmt.string_to_symbol("fcch")
+            fcch_tags = self.get_tags_in_range(0, start, stop, key)
+            if fcch_tags:
+                self.sch_offset = fcch_tags[0].offset + int(round(8*self.burst_size+0*self.guard_period)) #156.25 is number of GMSK symbols per timeslot, 
+                                                                                       #8.25 is arbitrary safety margin in order to avoid cutting boundary of SCH burst
+                self.state = self.states["reaching_sch_burst"]
+            
+        elif self.state == self.states["reaching_sch_burst"]:
+            samples_left = self.sch_offset-self.nitems_written(0)
+            if samples_left <= len(in0)-self.history():
+                to_consume = samples_left
+                self.state = self.states["sch_at_input_buffer"]
+
+        elif self.state == self.states["sch_at_input_buffer"]:
+            offset = self.nitems_written(0)
+            key = pmt.string_to_symbol("sch")
+            value =  pmt.from_double(0)
+            self.add_item_tag(0,offset, key, value)
+            self.state = self.states["waiting_for_fcch_tag"]
+            self.sch_receiver.get_chan_imp_resp(in0[0:self.block_size+self.guard_period])
+#            plot(unwrap(angle(in0[0:2*self.block_size])))
+#            show()
+
+        out[:] = in0[self.history()-1:]
+        return to_consume
+