Ok, I'll try to resume the tools fro breaking chaos modulation, starting with spicy

# Author: Travis Oliphant
	# 1999 -- 2002
	from __future__ import division, print_function, absolute_import
	import operator
	import threading
	import sys
	import timeit
	from . import sigtools, dlti
	from ._upfirdn import upfirdn, _output_len
	from scipy._lib.six import callable
	from scipy._lib._version import NumpyVersion
	from scipy import fftpack, linalg
	from numpy import (allclose, angle, arange, argsort, array, asarray,
	atleast_1d, atleast_2d, cast, dot, exp, expand_dims,
	iscomplexobj, mean, ndarray, newaxis, ones, pi,
	poly, polyadd, polyder, polydiv, polymul, polysub, polyval,
	product, r_, ravel, real_if_close, reshape,
	roots, sort, take, transpose, unique, where, zeros,
	zeros_like)
	import numpy as np
	import math
	from scipy.special import factorial
	from .windows import get_window
	from ._arraytools import axis_slice, axis_reverse, odd_ext, even_ext, const_ext
	from .filter_design import cheby1, _validate_sos
	from .fir_filter_design import firwin
	if sys.version_info.major >= 3 and sys.version_info.minor >= 5:
	from math import gcd
	else:
	from fractions import gcd
	__all__ = ['correlate', 'fftconvolve', 'convolve', 'convolve2d', 'correlate2d',
	'order_filter', 'medfilt', 'medfilt2d', 'wiener', 'lfilter',
	'lfiltic', 'sosfilt', 'deconvolve', 'hilbert', 'hilbert2',
	'cmplx_sort', 'unique_roots', 'invres', 'invresz', 'residue',
	'residuez', 'resample', 'resample_poly', 'detrend',
	'lfilter_zi', 'sosfilt_zi', 'sosfiltfilt', 'choose_conv_method',
	'filtfilt', 'decimate', 'vectorstrength']
	_modedict = {'valid': 0, 'same': 1, 'full': 2}
	_boundarydict = {'fill': 0, 'pad': 0, 'wrap': 2, 'circular': 2, 'symm': 1,
	'symmetric': 1, 'reflect': 4}
	_rfft_mt_safe = (NumpyVersion(np.__version__) >= '1.9.0.dev-e24486e')
	_rfft_lock = threading.Lock()
	def _valfrommode(mode):
	try:
	val = _modedict[mode]
	except KeyError:
	if mode not in [0, 1, 2]:
	raise ValueError("Acceptable mode flags are 'valid' (0),"
	" 'same' (1), or 'full' (2).")
	val = mode
	return val
	def _bvalfromboundary(boundary):
	try:
	val = _boundarydict[boundary] << 2
	except KeyError:
	if val not in [0, 1, 2]:
	raise ValueError("Acceptable boundary flags are 'fill', 'wrap'"
	" (or 'circular'), \n and 'symm'"
	" (or 'symmetric').")
	val = boundary << 2
	return val
	def _inputs_swap_needed(mode, shape1, shape2):
	"""
	If in 'valid' mode, returns whether or not the input arrays need to be
	swapped depending on whether `shape1` is at least as large as `shape2` in
	every dimension.
	This is important for some of the correlation and convolution
	implementations in this module, where the larger array input needs to come
	before the smaller array input when operating in this mode.
	Note that if the mode provided is not 'valid', False is immediately
	returned.
	"""
	if mode == 'valid':
	ok1, ok2 = True, True
	for d1, d2 in zip(shape1, shape2):
	if not d1 >= d2:
	ok1 = False
	if not d2 >= d1:
	ok2 = False
	if not (ok1 or ok2):
	raise ValueError("For 'valid' mode, one must be at least "
	"as large as the other in every dimension")
	return not ok1
	return False
	def correlate(in1, in2, mode='full', method='auto'):
	r"""
	Cross-correlate two N-dimensional arrays.
	Cross-correlate `in1` and `in2`, with the output size determined by the
	`mode` argument.
	Parameters
	----------
	in1 : array_like
	First input.
	in2 : array_like
	Second input. Should have the same number of dimensions as `in1`.
	mode : str {'full', 'valid', 'same'}, optional
	A string indicating the size of the output:
	``full``
	The output is the full discrete linear cross-correlation
	of the inputs. (Default)
	``valid``
	The output consists only of those elements that do not
	rely on the zero-padding. In 'valid' mode, either `in1` or `in2`
	must be at least as large as the other in every dimension.
	``same``
	The output is the same size as `in1`, centered
	with respect to the 'full' output.
	method : str {'auto', 'direct', 'fft'}, optional
	A string indicating which method to use to calculate the correlation.
	``direct``
	The correlation is determined directly from sums, the definition of
	correlation.
	``fft``
	The Fast Fourier Transform is used to perform the correlation more
	quickly (only available for numerical arrays.)
	``auto``
	Automatically chooses direct or Fourier method based on an estimate
	of which is faster (default). See `convolve` Notes for more detail.
	.. versionadded:: 0.19.0
	Returns
	-------
	correlate : array
	An N-dimensional array containing a subset of the discrete linear
	cross-correlation of `in1` with `in2`.
	See Also
	--------
	choose_conv_method : contains more documentation on `method`.
	Notes
	-----
	The correlation z of two d-dimensional arrays x and y is defined as::
	z[...,k,...] = sum[..., i_l, ...] x[..., i_l,...] * conj(y[..., i_l - k,...])
	This way, if x and y are 1-D arrays and ``z = correlate(x, y, 'full')`` then
	.. math::
	z[k] = (x * y)(k - N + 1)
	= \sum_{l=0}^{||x||-1}x_l y_{l-k+N-1}^{*}
	for :math:`k = 0, 1, ..., ||x|| + ||y|| - 2`
	where :math:`||x||` is the length of ``x``, :math:`N = \max(||x||,||y||)`,
	and :math:`y_m` is 0 when m is outside the range of y.
	``method='fft'`` only works for numerical arrays as it relies on
	`fftconvolve`. In certain cases (i.e., arrays of objects or when
	rounding integers can lose precision), ``method='direct'`` is always used.
	Examples
	--------
	Implement a matched filter using cross-correlation, to recover a signal
	that has passed through a noisy channel.
	>>> from scipy import signal
	>>> sig = np.repeat([0., 1., 1., 0., 1., 0., 0., 1.], 128)
	>>> sig_noise = sig + np.random.randn(len(sig))
	>>> corr = signal.correlate(sig_noise, np.ones(128), mode='same') / 128
	>>> import matplotlib.pyplot as plt
	>>> clock = np.arange(64, len(sig), 128)
	>>> fig, (ax_orig, ax_noise, ax_corr) = plt.subplots(3, 1, sharex=True)
	>>> ax_orig.plot(sig)
	>>> ax_orig.plot(clock, sig[clock], 'ro')
	>>> ax_orig.set_title('Original signal')
	>>> ax_noise.plot(sig_noise)
	>>> ax_noise.set_title('Signal with noise')
	>>> ax_corr.plot(corr)
	>>> ax_corr.plot(clock, corr[clock], 'ro')
	>>> ax_corr.axhline(0.5, ls=':')
	>>> ax_corr.set_title('Cross-correlated with rectangular pulse')
	>>> ax_orig.margins(0, 0.1)
	>>> fig.tight_layout()
	>>> fig.show()
	"""
	in1 = asarray(in1)
	in2 = asarray(in2)
	if in1.ndim == in2.ndim == 0:
	return in1 * in2
	elif in1.ndim != in2.ndim:
	raise ValueError("in1 and in2 should have the same dimensionality")
	# Don't use _valfrommode, since correlate should not accept numeric modes
	try:
	val = _modedict[mode]
	except KeyError:
	raise ValueError("Acceptable mode flags are 'valid',"
	" 'same', or 'full'.")
	# this either calls fftconvolve or this function with method=='direct'
	if method in ('fft', 'auto'):
	return convolve(in1, _reverse_and_conj(in2), mode, method)
	# fastpath to faster numpy.correlate for 1d inputs when possible
	if _np_conv_ok(in1, in2, mode):
	return np.correlate(in1, in2, mode)
	# _correlateND is far slower when in2.size > in1.size, so swap them
	# and then undo the effect afterward if mode == 'full'. Also, it fails
	# with 'valid' mode if in2 is larger than in1, so swap those, too.
	# Don't swap inputs for 'same' mode, since shape of in1 matters.
	swapped_inputs = ((mode == 'full') and (in2.size > in1.size) or
	_inputs_swap_needed(mode, in1.shape, in2.shape))
	if swapped_inputs:
	in1, in2 = in2, in1
	if mode == 'valid':
	ps = [i - j + 1 for i, j in zip(in1.shape, in2.shape)]
	out = np.empty(ps, in1.dtype)
	z = sigtools._correlateND(in1, in2, out, val)
	else:
	ps = [i + j - 1 for i, j in zip(in1.shape, in2.shape)]
	# zero pad input
	in1zpadded = np.zeros(ps, in1.dtype)
	sc = [slice(0, i) for i in in1.shape]
	in1zpadded[sc] = in1.copy()
	if mode == 'full':
	out = np.empty(ps, in1.dtype)
	elif mode == 'same':
	out = np.empty(in1.shape, in1.dtype)
	z = sigtools._correlateND(in1zpadded, in2, out, val)
	if swapped_inputs:
	# Reverse and conjugate to undo the effect of swapping inputs
	z = _reverse_and_conj(z)
	return z
	def _centered(arr, newshape):
	# Return the center newshape portion of the array.
	newshape = asarray(newshape)
	currshape = array(arr.shape)
	startind = (currshape - newshape) // 2
	endind = startind + newshape
	myslice = [slice(startind[k], endind[k]) for k in range(len(endind))]
	return arr[tuple(myslice)]
https://github.com/scipy/scipy/blob/master/scipy/signal/signaltools.py