# Source code for deltasigma._logsmooth

# -*- coding: utf-8 -*-
# _logsmooth.py
# Module providing the logsmooth function
# 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 logsmooth() function
"""

from __future__ import division

import numpy as np
from scipy.linalg import norm

from ._dbp import dbp

[docs]def logsmooth(X, inBin, nbin=8, n=3):
"""Smooth the FFT and convert it to dB.

**Parameters:**

X : (N,) ndarray
The FFT data.

inBin : int
The bin index of the input sine wave (if any).

nbin : int, optional
The number of bins on which the averaging will be performed,
used *before* 3*inBin

n : int, optional
Around the location of the input signal and its harmonics (up to the
third harmonic), don't average for n bins.

The logsmooth algorithm uses nbin bins from 0 to 3*inBin,
thereafter the bin sizes are increased by a factor 1.1,
staying less than 2^10.

For the :math:n sets of bins:
:math:inBin + i, 2*inBin + i ... n*inBin+i, where :math:i \\in [0,2]
don't do averaging. This way, the noise BW
and the scaling of the tone and its harmonics are unchanged.

.. note::

Unfortunately, harmonics above the nth appear smaller than they
really are because their energy is averaged over many bins.

**Returns:**

f, p : tuple of 1d-arrays
The bins and smoothed FFT, expressed in dB.

.. seealso::

* :func:plotSpectrum, convenience function to first call
:func:logsmooth and then plot on a logarithmic x-axis its return
values.

* :func:circ_smooth, smoothing algorithm suitable for linear
x-axis plotting.

.. plot::

import matplotlib.pyplot as plt
import numpy as np
from deltasigma import dbv, ds_hann, figureMagic, logsmooth
T = 2 #s
Fs = 231e3 #Hz
N = int(np.round(T*Fs, 0)) # FFT points
freq = .1e3
t = np.arange(N)/Fs
u0 = np.sin(2*np.pi*t*freq)
u0 = u0 + .01*u0**2+.001*u0**3+.005*u0**4
U = np.fft.fft(u0 * ds_hann(N))/(N/4)
f = np.linspace(0, Fs, N + 1)
f = f[:N/2 + 1]
plt.subplot(211)
plt.semilogx(f, dbv(U[:N/2 + 1]))
plt.hold(True)
inBin = np.round(freq/Fs*N)
fS, US = logsmooth(U, inBin)
plt.semilogx(fS*Fs, US, 'r', linewidth=2.5)
plt.xlim([f*Fs, Fs/2])
plt.ylabel('U(f) [dB]')
figureMagic(xRange=[100, 1e4], yRange=[-400, 0], name='Spectrum')
plt.subplot(212)
plt.loglog(fS[1:]*Fs, np.diff(fS*Fs))
plt.xlabel('f [Hz]')
plt.ylabel('Averaging interval [Hz]')
figureMagic(xRange=[100, 1e4])
plt.show()

"""
# preliminary sanitization of the input
if not np.prod(X.shape) == max(X.shape):
raise ValueError('Expected a (N,) or (N, 1)-shaped array.')
if len(X.shape) > 1:
X = np.squeeze(X)
inBin = int(inBin)

N = X.shape
N2 = int(np.floor(N/2))
f1 = int(inBin % nbin)
startbin = np.concatenate((np.arange(f1, inBin, nbin),
np.arange(inBin, inBin + 3)
))
i = 1 # my fix
while i < n: # n can be big and xrange is not in Python3
startbin = np.concatenate((startbin,
np.arange(startbin[-1] + 1,
(inBin + 1)*(i + 1) - 1, nbin),
(i + 1)*(inBin + 1) - 1 + np.arange(0, 3)
))
i = i + 1
# the following is my fix - REP?
startbin = np.concatenate((startbin, np.array((startbin[-1] + 1,))))
m = startbin[-1] + nbin
while m < N2 - 1:
startbin = np.concatenate((startbin, np.array((m,))))
nbin = np.min((nbin*1.1, 2**10))
m = int(np.round(m + nbin, 0))

stopbin = np.concatenate((startbin[1:] - 1, np.array((N2 - 1,))))
f = ((startbin + stopbin)/2)/N
p = np.zeros(f.shape)
for i in range(f.shape):
p[i] = dbp(norm(X[startbin[i]:stopbin[i] + 1])**2/
(stopbin[i] - startbin[i] + 1))
return f, p