Changed directory structure.
Corrected clock_offset_corrector (for some streange and yet unknown reason fractional resampler eats strem tags for some values of sps).
(this commit may contain some changes that are not described)
diff --git a/python/receiver/chirpz.py b/python/receiver/chirpz.py
new file mode 100644
index 0000000..3043c44
--- /dev/null
+++ b/python/receiver/chirpz.py
@@ -0,0 +1,488 @@
+# This program is public domain
+# Authors: Paul Kienzle, Nadav Horesh
+"""
+Chirp z-transform.
+
+CZT: callable (x,axis=-1)->array
+ define a chirp-z transform that can be applied to different signals
+ZoomFFT: callable (x,axis=-1)->array
+ define a Fourier transform on a range of frequencies
+ScaledFFT: callable (x,axis=-1)->array
+ define a limited frequency FFT
+
+czt: array
+ compute the chirp-z transform for a signal
+zoomfft: array
+ compute the Fourier transform on a range of frequencies
+scaledfft: array
+ compute a limited frequency FFT for a signal
+"""
+__all__ = ['czt', 'zoomfft', 'scaledfft']
+
+import math, cmath
+
+import numpy as np
+from numpy import pi, arange
+from scipy.fftpack import fft, ifft, fftshift
+
+class CZT:
+ """
+ Chirp-Z Transform.
+
+ Transform to compute the frequency response around a spiral.
+ Objects of this class are callables which can compute the
+ chirp-z transform on their inputs. This object precalculates
+ constants for the given transform.
+
+ If w does not lie on the unit circle, then the transform will be
+ around a spiral with exponentially increasing radius. Regardless,
+ angle will increase linearly.
+
+ The chirp-z transform can be faster than an equivalent fft with
+ zero padding. Try it with your own array sizes to see. It is
+ theoretically faster for large prime fourier transforms, but not
+ in practice.
+
+ The chirp-z transform is considerably less precise than the
+ equivalent zero-padded FFT, with differences on the order of 1e-11
+ from the direct transform rather than the on the order of 1e-15 as
+ seen with zero-padding.
+
+ See zoomfft for a friendlier interface to partial fft calculations.
+ """
+ def __init__(self, n, m=None, w=1, a=1):
+ """
+ Chirp-Z transform definition.
+
+ Parameters:
+ ----------
+ n: int
+ The size of the signal
+ m: int
+ The number of points desired. The default is the length of the input data.
+ a: complex
+ The starting point in the complex plane. The default is 1.
+ w: complex or float
+ If w is complex, it is the ratio between points in each step.
+ If w is float, it serves as a frequency scaling factor. for instance
+ when assigning w=0.5, the result FT will span half of frequncy range
+ (that fft would result) at half of the frequncy step size.
+
+ Returns:
+ --------
+ CZT:
+ callable object f(x,axis=-1) for computing the chirp-z transform on x
+ """
+ if m is None:
+ m = n
+ if w is None:
+ w = cmath.exp(-1j*pi/m)
+ elif type(w) in (float, int):
+ w = cmath.exp(-1j*pi/m * w)
+ else:
+ w = cmath.sqrt(w)
+ self.w, self.a = w, a
+ self.m, self.n = m, n
+
+ k = arange(max(m,n))
+ wk2 = w**(k**2)
+ nfft = 2**nextpow2(n+m-1)
+ self._Awk2 = (a**-k * wk2)[:n]
+ self._nfft = nfft
+ self._Fwk2 = fft(1/np.hstack((wk2[n-1:0:-1], wk2[:m])), nfft)
+ self._wk2 = wk2[:m]
+ self._yidx = slice(n-1, n+m-1)
+
+ def __call__(self, x, axis=-1):
+ """
+ Parameters:
+ ----------
+ x: array
+ The signal to transform.
+ axis: int
+ Array dimension to operate over. The default is the final
+ dimension.
+
+ Returns:
+ -------
+ An array of the same dimensions as x, but with the length of the
+ transformed axis set to m. Note that this is a view on a much
+ larger array. To save space, you may want to call it as
+ y = czt(x).copy()
+ """
+ x = np.asarray(x)
+ if x.shape[axis] != self.n:
+ raise ValueError("CZT defined for length %d, not %d" %
+ (self.n, x.shape[axis]))
+ # Calculate transpose coordinates, to allow operation on any given axis
+ trnsp = np.arange(x.ndim)
+ trnsp[[axis, -1]] = [-1, axis]
+ x = x.transpose(*trnsp)
+ y = ifft(self._Fwk2 * fft(x*self._Awk2, self._nfft))
+ y = y[..., self._yidx] * self._wk2
+ return y.transpose(*trnsp)
+
+
+def nextpow2(n):
+ """
+ Return the smallest power of two greater than or equal to n.
+ """
+ return int(math.ceil(math.log(n)/math.log(2)))
+
+
+def ZoomFFT(n, f1, f2=None, m=None, Fs=2):
+ """
+ Zoom FFT transform definition.
+
+ Computes the Fourier transform for a set of equally spaced
+ frequencies.
+
+ Parameters:
+ ----------
+ n: int
+ size of the signal
+ m: int
+ size of the output
+ f1, f2: float
+ start and end frequencies; if f2 is not specified, use 0 to f1
+ Fs: float
+ sampling frequency (default=2)
+
+ Returns:
+ -------
+ A CZT instance
+ A callable object f(x,axis=-1) for computing the zoom FFT on x.
+
+ Sampling frequency is 1/dt, the time step between samples in the
+ signal x. The unit circle corresponds to frequencies from 0 up
+ to the sampling frequency. The default sampling frequency of 2
+ means that f1,f2 values up to the Nyquist frequency are in the
+ range [0,1). For f1,f2 values expressed in radians, a sampling
+ frequency of 1/pi should be used.
+
+ To graph the magnitude of the resulting transform, use::
+
+ plot(linspace(f1,f2,m), abs(zoomfft(x,f1,f2,m))).
+
+ Use the zoomfft wrapper if you only need to compute one transform.
+ """
+ if m is None: m = n
+ if f2 is None: f1, f2 = 0., f1
+ w = cmath.exp(-2j * pi * (f2-f1) / ((m-1)*Fs))
+ a = cmath.exp(2j * pi * f1/Fs)
+ return CZT(n, m=m, w=w, a=a)
+
+def ScaledFFT(n, m=None, scale=1.0):
+ """
+ Scaled fft transform definition.
+
+ Similar to fft, where the frequency range is scaled by a factor 'scale' and
+ divided into 'm-1' equal steps. Like the FFT, frequencies are arranged
+ from 0 to scale*Fs/2-delta followed by -scale*Fs/2 to -delta, where delta
+ is the step size scale*Fs/m for sampling frequence Fs. The intended use is in
+ a convolution of two signals, each has its own sampling step.
+
+ This is equivalent to:
+
+ fftshift(zoomfft(x, -scale, scale*(m-2.)/m, m=m))
+
+ For example:
+
+ m,n = 10,len(x)
+ sf = ScaledFFT(n, m=m, scale=0.25)
+ X = fftshift(fft(x))
+ W = linspace(-8, 8*(n-2.)/n, n)
+ SX = fftshift(sf(x))
+ SW = linspace(-2, 2*(m-2.)/m, m)
+ plot(X,W,SX,SW)
+
+ Parameters:
+ ----------
+ n: int
+ Size of the signal
+ m: int
+ The size of the output.
+ Default: m=n
+ scale: float
+ Frequenct scaling factor.
+ Default: scale=1.0
+
+ Returns:
+ -------
+ function
+ A callable f(x,axis=-1) for computing the scaled FFT on x.
+ """
+ if m is None:
+ m = n
+ w = np.exp(-2j * pi / m * scale)
+ a = w**(m//2)
+ transform = CZT(n=n, m=m, a=a, w=w)
+ return lambda x, axis=-1: fftshift(transform(x, axis), axes=(axis,))
+
+def scaledfft(x, m=None, scale=1.0, axis=-1):
+ """
+ Partial with a frequency scaling.
+ See ScaledFFT doc for details
+
+ Parameters:
+ ----------
+ x: input array
+ m: int
+ The length of the output signal
+ scale: float
+ A frequency scaling factor
+ axis: int
+ The array dimension to operate over. The default is the
+ final dimension.
+
+ Returns:
+ -------
+ An array of the same rank of 'x', but with the size if
+ the 'axis' dimension set to 'm'
+ """
+ return ScaledFFT(x.shape[axis], m, scale)(x,axis)
+
+def czt(x, m=None, w=1.0, a=1, axis=-1):
+ """
+ Compute the frequency response around a spiral.
+
+ Parameters:
+ ----------
+ x: array
+ The set of data to transform.
+ m: int
+ The number of points desired. The default is the length of the input data.
+ a: complex
+ The starting point in the complex plane. The default is 1.
+ w: complex or float
+ If w is complex, it is the ratio between points in each step.
+ If w is float, it is the frequency step scale (relative to the
+ normal dft frquency step).
+ axis: int
+ Array dimension to operate over. The default is the final
+ dimension.
+
+ Returns:
+ -------
+ An array of the same dimensions as x, but with the length of the
+ transformed axis set to m. Note that this is a view on a much
+ larger array. To save space, you may want to call it as
+ y = ascontiguousarray(czt(x))
+
+ See zoomfft for a friendlier interface to partial fft calculations.
+
+ If the transform needs to be repeated, use CZT to construct a
+ specialized transform function which can be reused without
+ recomputing constants.
+ """
+ x = np.asarray(x)
+ transform = CZT(x.shape[axis], m=m, w=w, a=a)
+ return transform(x,axis=axis)
+
+def zoomfft(x, f1, f2=None, m=None, Fs=2, axis=-1):
+ """
+ Compute the Fourier transform of x for frequencies in [f1, f2].
+
+ Parameters:
+ ----------
+ m: int
+ The number of points to evaluate. The default is the length of x.
+ f1, f2: float
+ The frequency range. If f2 is not specified, the range 0-f1 is assumed.
+ Fs: float
+ The sampling frequency. With a sampling frequency of
+ 10kHz for example, the range f1 and f2 can be expressed in kHz.
+ The default sampling frequency is 2, so f1 and f2 should be
+ in the range 0,1 to keep the transform below the Nyquist
+ frequency.
+ x : array
+ The input signal.
+ axis: int
+ The array dimension the transform operates over. The default is the
+ final dimension.
+
+ Returns:
+ -------
+ array
+ The transformed signal. The fourier transform will be calculate
+ at the points f1, f1+df, f1+2df, ..., f2, where df=(f2-f1)/m.
+
+ zoomfft(x,0,2-2./len(x)) is equivalent to fft(x).
+
+ To graph the magnitude of the resulting transform, use::
+
+ plot(linspace(f1,f2,m), abs(zoomfit(x,f1,f2,m))).
+
+ If the transform needs to be repeated, use ZoomFFT to construct a
+ specialized transform function which can be reused without
+ recomputing constants.
+ """
+ x = np.asarray(x)
+ transform = ZoomFFT(x.shape[axis], f1, f2=f2, m=m, Fs=Fs)
+ return transform(x,axis=axis)
+
+
+def _test1(x,show=False,plots=[1,2,3,4]):
+ norm = np.linalg.norm
+
+ # Normal fft and zero-padded fft equivalent to 10x oversampling
+ over=10
+ w = np.linspace(0,2-2./len(x),len(x))
+ y = fft(x)
+ wover = np.linspace(0,2-2./(over*len(x)),over*len(x))
+ yover = fft(x,over*len(x))
+
+ # Check that zoomfft is the equivalent of fft
+ y1 = zoomfft(x,0,2-2./len(y))
+
+ # Check that zoomfft with oversampling is equivalent to zero padding
+ y2 = zoomfft(x,0,2-2./len(yover), m=len(yover))
+
+ # Check that zoomfft works on a subrange
+ f1,f2 = w[3],w[6]
+ y3 = zoomfft(x,f1,f2,m=3*over+1)
+ w3 = np.linspace(f1,f2,len(y3))
+ idx3 = slice(3*over,6*over+1)
+
+ if not show: plots = []
+ if plots != []:
+ import pylab
+ if 0 in plots:
+ pylab.figure(0)
+ pylab.plot(x)
+ pylab.ylabel('Intensity')
+ if 1 in plots:
+ pylab.figure(1)
+ pylab.subplot(311)
+ pylab.plot(w,abs(y),'o',w,abs(y1))
+ pylab.legend(['fft','zoom'])
+ pylab.ylabel('Magnitude')
+ pylab.title('FFT equivalent')
+ pylab.subplot(312)
+ pylab.plot(w,np.angle(y),'o',w,np.angle(y1))
+ pylab.legend(['fft','zoom'])
+ pylab.ylabel('Phase (radians)')
+ pylab.subplot(313)
+ pylab.plot(w,abs(y)-abs(y1)) #,w,np.angle(y)-np.angle(y1))
+ #pylab.legend(['magnitude','phase'])
+ pylab.ylabel('Residuals')
+ if 2 in plots:
+ pylab.figure(2)
+ pylab.subplot(211)
+ pylab.plot(w,abs(y),'o',wover,abs(y2),wover,abs(yover))
+ pylab.ylabel('Magnitude')
+ pylab.title('Oversampled FFT')
+ pylab.legend(['fft','zoom','pad'])
+ pylab.subplot(212)
+ pylab.plot(wover,abs(yover)-abs(y2),
+ w,abs(y)-abs(y2[0::over]),'o',
+ w,abs(y)-abs(yover[0::over]),'x')
+ pylab.legend(['pad-zoom','fft-zoom','fft-pad'])
+ pylab.ylabel('Residuals')
+ if 3 in plots:
+ pylab.figure(3)
+ ax1=pylab.subplot(211)
+ pylab.plot(w,abs(y),'o',w3,abs(y3),wover,abs(yover),
+ w[3:7],abs(y3[::over]),'x')
+ pylab.title('Zoomed FFT')
+ pylab.ylabel('Magnitude')
+ pylab.legend(['fft','zoom','pad'])
+ pylab.plot(w3,abs(y3),'x')
+ ax1.set_xlim(f1,f2)
+ ax2=pylab.subplot(212)
+ pylab.plot(wover[idx3],abs(yover[idx3])-abs(y3),
+ w[3:7],abs(y[3:7])-abs(y3[::over]),'o',
+ w[3:7],abs(y[3:7])-abs(yover[3*over:6*over+1:over]),'x')
+ pylab.legend(['pad-zoom','fft-zoom','fft-pad'])
+ ax2.set_xlim(f1,f2)
+ pylab.ylabel('Residuals')
+ if plots != []:
+ pylab.show()
+
+ err = norm(y-y1)/norm(y)
+ #print "direct err %g"%err
+ assert err < 1e-10, "error for direct transform is %g"%(err,)
+ err = norm(yover-y2)/norm(yover)
+ #print "over err %g"%err
+ assert err < 1e-10, "error for oversampling is %g"%(err,)
+ err = norm(yover[idx3]-y3)/norm(yover[idx3])
+ #print "range err %g"%err
+ assert err < 1e-10, "error for subrange is %g"%(err,)
+
+def _testscaled(x):
+ n = len(x)
+ norm = np.linalg.norm
+ assert norm(fft(x)-scaledfft(x)) < 1e-10
+ assert norm(fftshift(fft(x))[n/4:3*n/4] - fftshift(scaledfft(x,scale=0.5,m=n/2))) < 1e-10
+
+def test(demo=None,plots=[1,2,3]):
+ # 0: Gauss
+ t = np.linspace(-2,2,128)
+ x = np.exp(-t**2/0.01)
+ _test1(x, show=(demo==0), plots=plots)
+
+ # 1: Linear
+ x=[1,2,3,4,5,6,7]
+ _test1(x, show=(demo==1), plots=plots)
+
+ # Check near powers of two
+ _test1(range(126-31), show=False)
+ _test1(range(127-31), show=False)
+ _test1(range(128-31), show=False)
+ _test1(range(129-31), show=False)
+ _test1(range(130-31), show=False)
+
+ # Check transform on n-D array input
+ x = np.reshape(np.arange(3*2*28),(3,2,28))
+ y1 = zoomfft(x,0,2-2./28)
+ y2 = zoomfft(x[2,0,:],0,2-2./28)
+ err = np.linalg.norm(y2-y1[2,0])
+ assert err < 1e-15, "error for n-D array is %g"%(err,)
+
+ # 2: Random (not a test condition)
+ if demo==2:
+ x = np.random.rand(101)
+ _test1(x, show=True, plots=plots)
+
+ # 3: Spikes
+ t=np.linspace(0,1,128)
+ x=np.sin(2*pi*t*5)+np.sin(2*pi*t*13)
+ _test1(x, show=(demo==3), plots=plots)
+
+ # 4: Sines
+ x=np.zeros(100)
+ x[[1,5,21]]=1
+ _test1(x, show=(demo==4), plots=plots)
+
+ # 5: Sines plus complex component
+ x += 1j*np.linspace(0,0.5,x.shape[0])
+ _test1(x, show=(demo==5), plots=plots)
+
+ # 6: Scaled FFT on complex sines
+ x += 1j*np.linspace(0,0.5,x.shape[0])
+ if demo == 6:
+ demo_scaledfft(x,0.25,200)
+ _testscaled(x)
+
+
+def demo_scaledfft(v, scale, m):
+ import pylab
+ shift = pylab.fftshift
+ n = len(v)
+ x = pylab.linspace(-0.5, 0.5 - 1./n, n)
+ xz = pylab.linspace(-scale*0.5, scale*0.5*(m-2.)/m, m)
+ pylab.figure()
+ pylab.plot(x, shift(abs(fft(v))), label='fft')
+ pylab.plot(x, shift(abs(scaledfft(v))),'ro', label='x1 scaled fft')
+ pylab.plot(xz, abs(zoomfft(v, -scale, scale*(m-2.)/m, m=m)),
+ 'bo',label='zoomfft')
+ pylab.plot(xz, shift(abs(scaledfft(v, m=m, scale=scale))),
+ 'gx', label='x'+str(scale)+' scaled fft')
+ pylab.gca().set_yscale('log')
+ pylab.legend()
+ pylab.show()
+
+if __name__ == "__main__":
+ # Choose demo in [0,4] to show plot, or None for testing only
+ test(demo=None)
+
diff --git a/python/receiver/clock_offset_control.py b/python/receiver/clock_offset_control.py
new file mode 100644
index 0000000..506b0fd
--- /dev/null
+++ b/python/receiver/clock_offset_control.py
@@ -0,0 +1,102 @@
+#!/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 gnuradio import gr
+import pmt
+from threading import Timer
+
+class clock_offset_control(gr.basic_block):
+ """
+ docstring for block clock_offset_control
+ """
+ def __init__(self, fc, samp_rate):
+ gr.basic_block.__init__(self,
+ name="gsm_clock_offset_control",
+ in_sig=[],
+ out_sig=[])
+ self.fc = fc
+ self.samp_rate = samp_rate
+ self.message_port_register_in(pmt.intern("measurements"))
+ self.set_msg_handler(pmt.intern("measurements"), self.process_measurement)
+ self.message_port_register_out(pmt.intern("ppm"))
+ self.alfa = 0.3
+ self.ppm_estimate = -1e6
+ self.first_measurement = True
+ self.counter = 0
+ self.last_state = ""
+ self.timer = Timer(0.5, self.timed_reset)
+ self.last_ppm = -1e6
+
+ def process_measurement(self,msg):
+ if pmt.is_tuple(msg):
+ key = pmt.symbol_to_string(pmt.tuple_ref(msg,0))
+ if key == "freq_offset":
+ freq_offset = pmt.to_double(pmt.tuple_ref(msg,1))
+ ppm = -freq_offset/self.fc*1.0e6
+ state = pmt.symbol_to_string(pmt.tuple_ref(msg,2))
+
+ self.last_state = state
+
+ if abs(ppm) > 100: #safeguard against flawed measurements
+ ppm = 0
+ self.reset()
+
+ if state == "fcch_search":
+ msg_ppm = pmt.from_double(ppm)
+ self.message_port_pub(pmt.intern("ppm"), msg_ppm)
+ self.timer.cancel()
+ self.timer = Timer(0.5, self.timed_reset)
+ self.timer.start()
+ elif state == "synchronized":
+ self.timer.cancel()
+ if self.first_measurement:
+ self.ppm_estimate = ppm
+ self.first_measurement = False
+ else:
+ self.ppm_estimate = (1-self.alfa)*self.ppm_estimate+self.alfa*ppm
+
+ if self.counter == 5:
+ self.counter = 0
+ if abs(self.last_ppm-self.ppm_estimate) > 0.1:
+ msg_ppm = pmt.from_double(ppm)
+ self.message_port_pub(pmt.intern("ppm"), msg_ppm)
+ self.last_ppm = self.ppm_estimate
+ else:
+ self.counter=self.counter+1
+ elif state == "sync_loss":
+ self.reset()
+ msg_ppm = pmt.from_double(0.0)
+ self.message_port_pub(pmt.intern("ppm"), msg_ppm)
+
+
+ def timed_reset(self):
+ if self.last_state != "synchronized":
+# print "conditional reset"
+ self.reset()
+ msg_ppm = pmt.from_double(0.0)
+ self.message_port_pub(pmt.intern("ppm"), msg_ppm)
+
+ def reset(self):
+ self.ppm_estimate = -1e6
+ self.counter = 0
+ self.first_measurement = True
+
diff --git a/python/receiver/fcch_burst_tagger.py b/python/receiver/fcch_burst_tagger.py
new file mode 100644
index 0000000..56fead9
--- /dev/null
+++ b/python/receiver/fcch_burst_tagger.py
@@ -0,0 +1,130 @@
+#!/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 gsm.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.block_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.block_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
+
+ print "input_items:",len(in0)
+ print "coarse_idx",coarse_idx
+ print "coarse_idx+precise_idx",coarse_idx+precise_idx
+
+ 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(2*self.guard_period)):
+ 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/receiver/fcch_detector.py b/python/receiver/fcch_detector.py
new file mode 100644
index 0000000..627dd00
--- /dev/null
+++ b/python/receiver/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/receiver/receiver_hier.py b/python/receiver/receiver_hier.py
new file mode 100644
index 0000000..aa8fda3
--- /dev/null
+++ b/python/receiver/receiver_hier.py
@@ -0,0 +1,63 @@
+#!/usr/bin/env python
+
+import weakref
+import gsm
+from gnuradio.eng_option import eng_option
+from gnuradio import gr, gru, blocks
+from gnuradio import filter
+
+class receiver_hier(gr.hier_block2):
+ def __init__(self, input_rate, osr=4, arfcn=0):
+ gr.hier_block2.__init__(self,
+ "receiver_hier",
+ gr.io_signature(1, 1, gr.sizeof_gr_complex),
+ gr.io_signature(1, 1, 142*gr.sizeof_float))
+ #set rates
+ gsm_symb_rate = 1625000/6.0
+
+ self.message_port_register_hier_in("bursts")
+ self.message_port_register_hier_in("measurements")
+
+ self.input_rate = input_rate
+ self.osr = osr
+ self.arfcn = arfcn
+ self.sps = input_rate / (gsm_symb_rate * osr)
+
+ #create accompaning blocks
+ self.filtr = self._set_filter()
+ self.interpolator = self._set_interpolator()
+ self.receiver = self._set_receiver()
+ self.connect(self, self.filtr, self.interpolator, self.receiver, self)
+# self.connect(self, self.interpolator, self.receiver, self)
+ self.msg_connect(self.receiver, "bursts", weakref.proxy(self), "bursts")
+ self.msg_connect(self.receiver, "measurements", weakref.proxy(self), "measurements")
+
+ def _set_filter(self):
+ filter_cutoff = 125e3
+ filter_t_width = 10e3
+ offset = 0
+
+ filter_taps = filter.firdes.low_pass(1.0, self.input_rate, filter_cutoff, filter_t_width, filter.firdes.WIN_HAMMING)
+ filtr = filter.freq_xlating_fir_filter_ccf(1, filter_taps, offset, self.input_rate)
+ return filtr
+
+ def _set_interpolator(self):
+ interpolator = filter.fractional_resampler_cc(0, self.sps)
+ return interpolator
+
+ def _set_receiver(self):
+ receiver = gsm.receiver(self.osr, self.arfcn)
+ return receiver
+
+ def set_center_frequency(self, center_freq):
+ self.filtr.set_center_freq(center_freq)
+
+ def set_timing(self, timing_offset):
+ pass
+
+ def set_arfcn(self,arfcn):
+ self.receiver.set_arfcn(arfcn)
+
+ def reset(self):
+ self.receiver.reset()
+
diff --git a/python/receiver/sch_detector.py b/python/receiver/sch_detector.py
new file mode 100644
index 0000000..ca3c423
--- /dev/null
+++ b/python/receiver/sch_detector.py
@@ -0,0 +1,156 @@
+#!/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.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,self.L+1,mode='constant',axis=0)
+
+ print "correlation_bl.argmax()",argmax(abs(correlation_bl))
+ print "power_bl_mov_avg.argmax()",(power_bl_mov_avg).argmax()
+ print 'unravel_index(correlation_bl.argmax(), correlation_bl.shape)',unravel_index(argmax(abs(correlation_bl)), 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)
+ (r_corrmax, c_corrmax)=unravel_index(argmax(abs(correlation_bl)), correlation_bl.shape)
+ (r_powmax, c_powmax)=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):
+# if ii == c_powmax:
+# plot(range(ii,len(correlation_bl[:,0])*self.OSR,self.OSR),power_bl_mov_avg[:,ii]*5e6,'g.')
+# else:
+# plot(range(ii,len(correlation_bl[:,0])*self.OSR,self.OSR),power_bl_mov_avg[:,ii]*5e6,'r.')
+# show()
+# figure()
+ print 'r_powmax: ',r_powmax
+# plot(abs(correlation_bl[range(r_powmax-(self.L+1)/2+1,r_powmax+(self.L+1)/2+1), c_powmax]),'g')
+# hold(True)
+# plot(abs(correlation_bl[range(r_corrmax-(self.L+1)/2+1,r_corrmax+(self.L+1)/2+1), c_corrmax]),'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
+
+ def get_OSR(self):
+ return self.OSR
+
+ def set_OSR(self, OSR):
+ self.OSR = OSR
+ 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.sch_receiver = sch_receiver(OSR)
+