# Source code for deltasigma._peakSNR

# -*- coding: utf-8 -*-
# _peakSNR.py
# Module providing peakSNR
# This file is part of python-deltasigma.
#
# python-deltasigma is a 1:1 Python replacement of Richard Schreier's
# MATLAB delta sigma toolbox (aka "delsigma"), upon which it is heavily based.
# The delta sigma toolbox is (c) 2009, Richard Schreier.
#
# python-deltasigma 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
# LICENSE file for the licensing terms.

"""Module providing the peakSNR function.
"""

from __future__ import division

import numpy as np
from numpy.linalg import lstsq

from ._config import _debug
from ._dbv import dbv

[docs]def peakSNR(snr, amp):
"""Find the SNR peak by fitting the SNR curve.

A smooth curve is fitted to the top end of the SNR vs input amplitude
data with Legendre's least-squares method.

The curve fitted to the data is:

.. math::

y = ax + \\frac{b}{x - c}

**Parameters:**

snr : 1D ndarray or array-like
Signal to Noise Ratio array, expressed in decibels (dB).

amp : 1D ndarray or array-like
Amplitude array, expressed in decibels (dB).

The two arrays snr and amp should have the same size.

**Returns:**

(peak_snr, peak_amp) : tuple
The peak SNR value and its corresponding amplitude,
expressed in dB.
"""

# sanitize inputs
snr = np.atleast_1d(np.asarray(snr))
amp = np.atleast_1d(np.asarray(amp))
# Delete garbage data
for check in (np.isinf, np.isnan):
i = check(snr)
snr = np.delete(snr, np.where(i))
amp = np.delete(amp, np.where(i))
for x in (snr, amp):
if len(x.shape) > 1 and np.prod(x.shape) != max(x.shape):
raise ValueError("snr and amp should be vectors (ndim=1)" +
"snr.shape: %s, amp.shape: %s" % (snr.shape, amp.shape))
# All good
max_snr = np.max(snr)
i = np.flatnonzero(snr > max_snr - 10)
min_i = np.min(i)
max_i = np.max(i)
j = np.flatnonzero(snr[min_i:max_i + 1] < max_snr - 15)
if j:
max_i = min_i + np.min(j) - 2
i = np.arange(min_i, max_i + 1)
snr = 10.0**(snr[i]/20)
amp = 10.0**(amp[i]/20)
#n = max(i.shape) # unused variable, REP?
c = np.max(amp)*1.05
# fit y = ax + b/(x-c) to the data
A = np.hstack((amp.reshape(-1, 1), 1.0/(amp.reshape(-1, 1) - c)))
ab = np.linalg.lstsq(A, snr.reshape((-1, 1)))
peak_amp = c - np.sqrt(ab[1, 0]/ab[0, 0])
peak_snr = np.dot(np.array([[peak_amp, 1./(peak_amp-c)]]), ab) #None check mcode
peak_snr = dbv(peak_snr)
peak_amp = dbv(peak_amp)
if _debug:
import pylab as plt
pred = np.dot(A, ab)
hold = plt.ishold()
plt.hold(True)
plt.plot(dbv(amp), dbv(pred), '-', color='b')
plt.hold(hold)
return peak_snr, peak_amp