# -*- coding: utf-8 -*-
#
# Modulation Error Rate
#
# http://www.opendigitalradio.org
# Licence: The MIT License, see notice at the end of this file
import datetime
import os
import logging
import time
try:
logging_path = os.path.dirname(logging.getLoggerClass().root.handlers[0].baseFilename)
except:
logging_path = "/tmp/"
import numpy as np
import src.const
import scipy
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt
class Symbol_align:
"""
Find the phase offset to the start of the DAB symbols in an
unaligned dab signal.
"""
def __init__(self, sample_rate):
self.c = src.const.const(sample_rate)
pass
def _calc_offset_to_first_symbol_without_prefix(self, tx, debug=False):
tx_orig = tx[0:-self.c.T_U]
tx_cut_prefix = tx[self.c.T_U:]
tx_product = np.abs(tx_orig - tx_cut_prefix)
tx_product_avg = np.correlate(
tx_product,
np.ones(self.c.T_C),
mode='valid')
tx_product_avg_min_filt = \
scipy.ndimage.filters.minimum_filter1d(
tx_product_avg,
int(1.5 * self.c.T_S)
)
peaks = np.ravel(np.where(tx_product_avg == tx_product_avg_min_filt))
offset = peaks[np.argmin([tx_product_avg[peak] for peak in peaks])]
if debug:
fig = plt.figure(figsize=(9, 9))
ax = fig.add_subplot(4, 1, 1)
ax.plot(tx_product)
ylim = ax.get_ylim()
for peak in peaks:
ax.plot((peak, peak), (ylim[0], ylim[1]))
if peak == offset:
ax.text(peak, ylim[0] + 0.3 * np.diff(ylim), "offset", rotation=90)
else:
ax.text(peak, ylim[0] + 0.2 * np.diff(ylim), "peak", rotation=90)
ax.set_xlabel("Sample")
ax.set_ylabel("Conj. Product")
ax.set_title("Difference with shifted self")
ax = fig.add_subplot(4, 1, 2)
ax.plot(tx_product_avg)
ylim = ax.get_ylim()
for peak in peaks:
ax.plot((peak, peak), (ylim[0], ylim[1]))
if peak == offset:
ax.text(peak, ylim[0] + 0.3 * np.diff(ylim), "offset", rotation=90)
else:
ax.text(peak, ylim[0] + 0.2 * np.diff(ylim), "peak", rotation=90)
ax.set_xlabel("Sample")
ax.set_ylabel("Conj. Product")
ax.set_title("Moving Average")
ax = fig.add_subplot(4, 1, 3)
ax.plot(tx_product_avg_min_filt)
ylim = ax.get_ylim()
for peak in peaks:
ax.plot((peak, peak), (ylim[0], ylim[1]))
if peak == offset:
ax.text(peak, ylim[0] + 0.3 * np.diff(ylim), "offset", rotation=90)
else:
ax.text(peak, ylim[0] + 0.2 * np.diff(ylim), "peak", rotation=90)
ax.set_xlabel("Sample")
ax.set_ylabel("Conj. Product")
ax.set_title("Min Filter")
ax = fig.add_subplot(4, 1, 4)
tx_product_crop = tx_product[peaks[0]-50:peaks[0]+50]
x = range(tx_product.shape[0])[peaks[0]-50:peaks[0]+50]
ax.plot(x, tx_product_crop)
ylim = ax.get_ylim()
ax.plot((peaks[0], peaks[0]), (ylim[0], ylim[1]))
ax.set_xlabel("Sample")
ax.set_ylabel("Conj. Product")
ax.set_title("Difference with shifted self")
fig.tight_layout()
# "offset" measures where the shifted signal matches the
# original signal. Therefore we have to subtract the size
# of the shift to find the offset of the symbol start.
return (offset + self.c.T_C) % self.c.T_S
def _remove_outliers(self, x, stds=5):
deviation_from_mean = np.abs(x - np.mean(x))
inlier_idxs = deviation_from_mean < stds * np.std(x)
x = x[inlier_idxs]
return x
def _calc_delta_angle(self, fft, debug=False):
# Introduce invariance against carrier
angles = np.angle(fft) % (np.pi / 2.)
# Calculate Angle difference and compensate jumps
deltas_angle = np.diff(angles)
deltas_angle[deltas_angle > np.pi/4.] =\
deltas_angle[deltas_angle > np.pi/4.] - np.pi/2.
deltas_angle[-deltas_angle > np.pi/4.] = \
deltas_angle[-deltas_angle > np.pi/4.] + np.pi/2.
deltas_angle = self._remove_outliers(deltas_angle)
delta_angle = np.mean(deltas_angle)
delta_angle_std = np.std(deltas_angle)
if debug:
plt.subplot(211)
plt.plot(angles, 'p')
plt.subplot(212)
plt.plot(deltas_angle, 'p')
return delta_angle
def _delta_angle_to_samples(self, angle):
return - angle / self.c.phase_offset_per_sample
def _calc_sample_offset(self, sig, debug=False):
assert sig.shape[0] == self.c.T_U,\
"Input length is not a Symbol without cyclic prefix"
fft = np.fft.fftshift(np.fft.fft(sig))
fft_crop = np.delete(fft[self.c.FFT_start:self.c.FFT_end], self.c.FFT_delete)
delta_angle = self._calc_delta_angle(fft_crop, debug=debug)
delta_sample = self._delta_angle_to_samples(delta_angle)
delta_sample_int = np.round(delta_sample).astype(int)
error = np.abs(delta_sample_int - delta_sample)
if error > 0.1:
raise RuntimeError("Could not calculate sample offset")
return delta_sample_int
def calc_offset(self, tx):
off_sym = self._calc_offset_to_first_symbol_without_prefix(
tx, debug=False)
off_sam = self._calc_sample_offset(
tx[off_sym:off_sym + self.c.T_U])
off = (off_sym + off_sam) % self.c.T_S
assert self._calc_sample_offset(tx[off:off + self.c.T_U]) == 0, \
"Failed to calculate offset"
return off
def crop_symbol_without_cyclic_prefix(self, tx):
off = self.calc_offset(tx)
return tx[
off:
off+self.c.T_U
]
# The MIT License (MIT)
#
# Copyright (c) 2017 Andreas Steger
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.