sigp.pyold

import math
import warnings

import numpy as np
from numpy import ma
import scipy as sp
from numpy import linalg as la 
import matplotlib
rcParams = matplotlib.rcParams

import matplotlib.cbook as cbook
from matplotlib.cbook import _string_to_bool, mplDeprecation
import matplotlib.collections as mcoll
import matplotlib.colors as mcolors
import matplotlib.contour as mcontour
import matplotlib.dates as _  # <-registers a date unit converter
from matplotlib import docstring
import matplotlib.image as mimage
import matplotlib.legend as mlegend
import matplotlib.lines as mlines
import matplotlib.markers as mmarkers
import matplotlib.mlab as mlab
import matplotlib.path as mpath
import matplotlib.patches as mpatches
import matplotlib.quiver as mquiver
import matplotlib.stackplot as mstack
import matplotlib.streamplot as mstream
import matplotlib.table as mtable
import matplotlib.text as mtext
import matplotlib.ticker as mticker
import matplotlib.transforms as mtransforms
import matplotlib.tri as mtri
import matplotlib.transforms as mtrans
from matplotlib.container import BarContainer, ErrorbarContainer, StemContainer
from matplotlib.axes._base import _AxesBase
from matplotlib.axes._base import _process_plot_format
import pylab as pl

def hanning(x):
	"""
	w=hanning(n)
	Return the n point hanning window
	x_windows=hanning(x)
	Returns x as the hanning windowing array x_windowed
	The windowed signal is then x*x_window
	
	"""
	if isinstance(x,(list, tuple, np.ndarray)):
		n=np.max(x.shape)
		f=hanning(n)
		f,dud=np.meshgrid(f[0,:],np.arange(x.shape[0]))
	else:

		n=x
		f=np.reshape((np.sin(np.pi*np.arange(n)/(n-1))**2)*np.sqrt(8/3),(1,-1))
		f=f/np.linalg.norm(f)*np.sqrt(n)


	return f
def blackwin(x):
	"""
	w=blackwin(n)
	Return the n point Blackman window
	x_windows=blackwin(x)
	Returns x as the Blackman windowing array x_window
	The windowed signal is then x*x_window
	"""
	disp('untested')
	if isinstance(x,(list, tuple, np.ndarray)):
		n=np.max(x.shape)
		f=hanning(n)
		f,dud=np.meshgrid(f[0,:],np.arange(x.shape[0]))
	else:

		n=x
		f=np.reshape((0.42-0.5*np.cos(2*np.pi*(np.arange(n)+.5))/(n)+.08*np.cos(4*np.pi*(np.arange(n)+.5))/(n))*np.sqrt(5000/1523),(1,-1))
		f=f/np.linalg.norm(f)*np.sqrt(n)


	return f
def expwin(x,ts=.75):
	"""
	w=expwin(n)
	Return the n point exponential window
	x_windows=expwin(x)
	Returns x as the expwin windowing array x_windowed
	The windowed signal is then x*x_window
	The optional second argument set the 5% "settling time" of the window. Default is ts=0.75 
	"""
	print('untested')
	tc=-ts/np.log(.05)
	if isinstance(x,(list, tuple, np.ndarray)):
		n=np.max(x.shape)
		f=hanning(n)
		f,dud=np.meshgrid(f[0,:],np.arange(x.shape[0]))
	else:
		n=x
		v=(n-1)/n*np.arange(n)+(n-1)/n/2
		f=exp(-v/tc/(n-1))
		f=f/np.linalg.norm(f)*np.sqrt(n)    
		f=np.reshape(f,(1,-1))
		f=f/np.linalg.norm(f)*np.sqrt(n)


	return f
def hammwin(x):
	"""
	w=hammwin(n)
	Return the n point hamming window
	x_windows=hamming(x)
	Returns x as the hamming windowingarray x_windowed
	The windowed signal is then x*x_window
	"""
	print('untested')
	if isinstance(x,(list, tuple, np.ndarray)):
		n=np.max(x.shape)
		f=hanning(n)
		f,dud=np.meshgrid(f[0,:],np.arange(x.shape[0]))
	else:

		n=x
		f=np.reshape((0.54-0.46*np.cos(2*np.pi*(np.arange(n))/(n-1)))*np.sqrt(5000/1987),(1,-1))
		f=f/np.linalg.norm(f)*np.sqrt(n)


	return f
def flatwin(x):
	"""
	w=flatwin(n)
	Return the n point flat top window
	x_windows=flatwin(x)
	Returns x as the flat top windowing array x_windowed
	The windowed signal is then x*x_window
	McConnell, K. G., "Vibration Testing: Theory and Practice," Wiley, 1995.
	"""
	print('untested')
	if isinstance(x,(list, tuple, np.ndarray)):
		n=np.max(x.shape)
		f=hanning(n)
		f,dud=np.meshgrid(f[0,:],np.arange(x.shape[0]))
	else:

		n=x
		f=np.reshape((1.0-1.933*np.cos(2*np.pi*(np.arange(n))/(n-1))+1.286*np.cos(4*np.pi*(np.arange(n))/(n-1))-0.338*np.cos(6*np.pi*(np.arange(n))/(n-1))+0.032*np.cos(8*np.pi*(np.arange(n))/(n-1))),(1,-1))
		f=f/np.linalg.norm(f)*np.sqrt(n)
		
	return f
def boxwin(x):
	"""
	w=boxwin(n)
	Return the n point box window (uniform)
	x_windows=boxwin(x)
	Returns x as the boxwin windowing array x_windowed
	The windowed signal is then x*x_window
	"""
	print('untested')
	if isinstance(x,(list, tuple, np.ndarray)):
		n=np.max(x.shape)
		f=hanning(n)
		f,dud=np.meshgrid(f[0,:],np.arange(x.shape[0]))
	else:

		n=x
		#f=np.reshape((1.0-1.933*np.cos(2*np.pi*(np.arange(n))/(n-1))+1.286*np.cos(4*np.pi*(np.arange(n))/(n-1))-0.338*np.cos(6*np.pi*(np.arange(n))/(n-1))+0.032*np.cos(8*np.pi*(np.arange(n))/(n-1))),(1,-1))
		f=np.reshape( np.ones((1,n)),(1,-1))
		f=f/np.linalg.norm(f)*np.sqrt(n)
		
	return f


	
def hannwin(x):
	f=hanning(x)
	return f
def asd(x,t,window="hanning",ave=bool(True)):
	"""
	freq, Pxx=asd(x,t)
	Calculate the autospectrum (power spectrum) density of a signal x

	Examples
    --------
    >>> from scipy import signal
	>>> import numpy as np
    >>> import matplotlib.pyplot as plt
	>>> import vttools at vt

    Generate a 5 second test signal, a 10 V sine wave at 50 Hz, corrupted by
    0.001 V**2/Hz of white noise sampled at 1 kHz.

    >>> sample_freq = 1e3
    >>> tfinal = 5
	>>> fs=1000
	>>> A=10
    >>> freq = 50
    >>> noise_power = 0.001 * sample_freq / 2
    >>> time = np.arange(0,tfinal,1/fsample)
    >>> x = A*np.sin(2*np.pi*freq*time)
    >>> x += np.random.normal(scale=np.sqrt(noise_power),size=len(time)) # Adds noise to x

    Compute and plot the autopectrum density.

    >>> freq_vec, Pxx = vt.asd(x,time)
    >>> plt.semilogy(freq_vec, 20*np.log10(Pxx))
    >>> #plt.ylim([1e-7, 1e2])
    >>> plt.xlabel('frequency [Hz]')
    >>> plt.ylabel('PSD [V**2/Hz]')
    >>> plt.show()

    If we average the last half of the spectral density, to exclude the
    peak, we can recover the noise power on the signal.

    >>> np.mean(Pxx_den[256:])
    0.0009924865443739191

    Now compute and plot the power spectrum.

    >>> f, Pxx_spec = signal.periodogram(x, fs, 'flattop', scaling='spectrum')
    >>> plt.figure()
    >>> plt.semilogy(f, np.sqrt(Pxx_spec))
    >>> plt.ylim([1e-4, 1e1])
    >>> plt.xlabel('frequency [Hz]')
    >>> plt.ylabel('Linear spectrum [V RMS]')
    >>> plt.show()
	"""
	f, Pxx=crsd(x,x,t,window,ave)
	Pxx=Pxx.real
	return f, Pxx

def crsd(x,y,t,window="hanning",ave=bool(True)):
    print(t.shape)
	if t.shape[1]>1:
		t=t[0,2]-t[0,1]
	elif t.shape[0]>1:
		t=t[2,0]-t[1,0]
		disp('t must be a scalar or size (1,n)')
	  
	n=x.shape[1];
	
	#print(n)
	#print('y shape')
	#print(y.shape)
	import matplotlib.pyplot as plt
	if window=="none":
		a=1
	else:
		win=1
		if window=="hanning":#BLACKWIN, BOXWIN, EXPWIN, HAMMWIN, and TRIWIN
			win=hanning(x)
		elif window=="blackwin":
			win=blackwin(x)
		elif window=="boxwin":
			win=boxwin(x)
		elif window=="expwin":
			win=expwin(x)
		elif window=="hammwin":
			win=hamming(x)
		elif window=="triwin":
			win=triwin(x)
		elif window=="flatwin":
			win=triwin(x)

		y=y*win
		x=x*win
	
	ffty=np.fft.rfft(y,n,1)*t
	
	fftx=np.fft.rfft(x,n,1)*t

	Pxy=ffty*np.conj(fftx)/(n*t)*2
	#print('pxy shape')
	#print(Pxy.shape)

	if Pxy.shape[1]>1 and ave:
		Pxy=np.mean(Pxy,0)
		Pxy=np.reshape(Pxy,(1,-1))
	
	nfreq=1/t/2;
	#print(Pxy.shape)
	#print(Pxy.shape[1])
	f=np.reshape(np.linspace(0, nfreq, Pxy.shape[1]),(1,-1))

	return f, Pxy
def frfest(x,f,dt):#,n,options)
	"""
	f,Txf = tfest(x,f,dt) estimates the 
	H1 Frequency Response Function (FRF) between X and F.
	TFEST calculates the H1 FRF (Sxf/Sff).
	Freq is the frequency vector in Hertz and Txf is the
	Transfer Function in complex form.
	dt is the time step of the sampled data. 
	If dt is replaced by the time vector, dt will be extracted 
	from the time vector using dt = t[2]-t[1].

	NOT ENABLED YET
	If X and 
	F are matrices, TFEST  will find the frequency response
	function for each column and average the results unless 
	OPTIONS{1} is set to 'no'. 

	NOT ENABLED YET
	TFEST(X,Y,DT,N,OPTIONS) plots the Frequency Response Function 
	if there are no output arguments.  Click in the region of 
	interest to zoom in.  Each click will double the size of the 
	plot.  Double click to return to full scale.

	See also COH, ASD, CRSD, and TFPLOT.

	1994 by Joseph C. Slater
	Modifications:
	--------------
	7/6/00: Changed default FRF calculation from H2 to H1
		 Added H1, H2, and Hv options.
	4/13/15: Converted to Python
	"""

	if len(dt)>1:
		dt=dt(2)-dt(1)

	#if frftype==1
	freq,Pff=asd(f,dt);
	freq,Pxf=crsd(x,f,dt);
	tfunc=Pxf/Pff;
	## elseif frftype==2
	##   [freq,Pxx]=asd(x,dt,n,'noave');
	##   [freq,Pxf]=crsd(x,f,dt,n,'noave');
	##   tfunc=Pxx./Pxf;
	## elseif frftype==3
	##   [freq,Pxx]=asd(x,dt,n,'noave');
	##   [freq,Pff]=asd(f,dt,n,'noave');
	##   [freq,Pxf]=crsd(x,f,dt,n,'noave');
	##   for i=1:size(Pxx,1)
	##     for j=1:size(Pxx,2)
	##       frfm=[Pff(i,j) Pxf(i,j);Pxf(i,j) Pxx(i,j)];
	##       [v,d]=eig(frfm);
	##       [y,yi]=sort(diag(d));
	##       tfunc(i,j)=-v(yi(1),1)/v(yi(1),2);		
	##     end
	##   end
	## end

	tfunc=np.conj(tfunc)

	## sfft=size(Pxf);
	## if sfft(2)~=1 & ~strcmp(ave,'no')
	##   tfunc=mean(tfunc')';
	## end

	## % If no left hand arguments then plot results
	## if nargout==0
	##   subplot(211)
	##   plot(freq,20*log10(abs(tfunc)))

	##   xlabel('Frequency (Hz)')
	##   ylabel('Mag (dB)')
	##   grid
	##   zoom on
	##   subplot(212)
	##   phase=unwrap(angle(tfunc))*180/pi;
	##   %phase=angle(tfunc)*180/pi;
	##   plot(freq,phase)

	##   xlabel('Frequency (Hz)')
	##   ylabel('Phase (deg)')
	##   grid
	##   phmin_max=[floor(min(phase)/45)*45 ceil(max(phase)/45)*45];
	##   set(gca,'YLim',phmin_max)
	##   gridmin_max=round(phmin_max/90)*90;
	##   set(gca,'YTick',gridmin_max(1):90:gridmin_max(2))
	##   % set(gca,'GridLineStyle','--')
	##   % gridmin_max=round(phmin_max/45)*45;
	##   % set(gca,'YTick',gridmin_max(1):45:gridmin_max(2))
	##   set(gca,'GridLineStyle',':')
	##   set(gca,'YTickLabel',gridmin_max(1):90:gridmin_max(2))
	##   zoom on
	##   return
	## end

	## freqout=freq;
	## tfout=tfunc;
	return freq,tfunc
def frfestH1(x,f,dt):#,n,options)
	"""
	f,Txf = tfestH1(x,f,dt) estimates the 
	H1 Frequency Response Function (FRF) between x andf.
	tfestH1 calculates the H1 FRF (Sxf/Sff).
	f is the frequency vector in Hertz and Txf is the
	Transfer Function in complex form. 
	If dt is replaced by the time vector, dt will be extracted 
	from the time vector using dt = t[2] - t[1].

	NOT ENABLED YET
	If X and 
	F are matrices, TFEST  will find the frequency response
	function for each column and average the results unless 
	OPTIONS{1} is set to 'no'. 

	NOT ENABLED YET
	TFEST(X,Y,DT,N,OPTIONS) plots the Frequency Response Function 
	if there are no output arguments.  Click in the region of 
	interest to zoom in.  Each click will double the size of the 
	plot.  Double click to return to full scale.

	See also COH, ASD, CRSD, and TFPLOT.

	1994 by Joseph C. Slater
	Modifications:
	--------------
	7/6/00: Changed default FRF calculation from H2 to H1
		 Added H1, H2, and Hv options.
	4/13/15: Converted to Python
	"""

	#if frftype==1
	print(f)
	print(dt)
	freq,Pff=asd(f,dt);
	freq,Pxf=crsd(x,f,dt);
	tfunc=Pxf/Pff;
	## elseif frftype==2
	##   [freq,Pxx]=asd(x,dt,n,'noave');
	##   [freq,Pxf]=crsd(x,f,dt,n,'noave');
	##   tfunc=Pxx./Pxf;
	## elseif frftype==3
	##   [freq,Pxx]=asd(x,dt,n,'noave');
	##   [freq,Pff]=asd(f,dt,n,'noave');
	##   [freq,Pxf]=crsd(x,f,dt,n,'noave');
	##   for i=1:size(Pxx,1)
	##     for j=1:size(Pxx,2)
	##       frfm=[Pff(i,j) Pxf(i,j);Pxf(i,j) Pxx(i,j)];
	##       [v,d]=eig(frfm);
	##       [y,yi]=sort(diag(d));
	##       tfunc(i,j)=-v(yi(1),1)/v(yi(1),2);		
	##     end
	##   end
	## end

	tfunc=np.conj(tfunc)

	## sfft=size(Pxf);
	## if sfft(2)~=1 & ~strcmp(ave,'no')
	##   tfunc=mean(tfunc')';
	## end

	## % If no left hand arguments then plot results
	## if nargout==0
	##   subplot(211)
	##   plot(freq,20*log10(abs(tfunc)))

	##   xlabel('Frequency (Hz)')
	##   ylabel('Mag (dB)')
	##   grid
	##   zoom on
	##   subplot(212)
	##   phase=unwrap(angle(tfunc))*180/pi;
	##   %phase=angle(tfunc)*180/pi;
	##   plot(freq,phase)

	##   xlabel('Frequency (Hz)')
	##   ylabel('Phase (deg)')
	##   grid
	##   phmin_max=[floor(min(phase)/45)*45 ceil(max(phase)/45)*45];
	##   set(gca,'YLim',phmin_max)
	##   gridmin_max=round(phmin_max/90)*90;
	##   set(gca,'YTick',gridmin_max(1):90:gridmin_max(2))
	##   % set(gca,'GridLineStyle','--')
	##   % gridmin_max=round(phmin_max/45)*45;
	##   % set(gca,'YTick',gridmin_max(1):45:gridmin_max(2))
	##   set(gca,'GridLineStyle',':')
	##   set(gca,'YTickLabel',gridmin_max(1):90:gridmin_max(2))
	##   zoom on
	##   return
	## end

	## freqout=freq;
	## tfout=tfunc;
	return freq,tfunc
def frfestH2(x,f,dt):#,n,options)
	"""
	f,Txf = tfestH2(X,F,DT,N,OPTIONS) estimates the 
	H2 Frequency Response Function (FRF) between X and F.
	tfestH2 calculates the 21 FRF (Sxx/Sxf).
	Freq is the frequency vector in Hertz and Txf is the
	Transfer Function in complex form.
	dt is the time step of the sampled data. 
	If dt is replaced by the time vector, dt will be extracted 
	from the time vector using dt = t[2]-t[1].

	NOT ENABLED YET
	If X and 
	F are matrices, TFEST  will find the frequency response
	function for each column and average the results unless 
	OPTIONS{1} is set to 'no'. 

	NOT ENABLED YET
	TFEST(X,Y,DT,N,OPTIONS) plots the Frequency Response Function 
	if there are no output arguments.  Click in the region of 
	interest to zoom in.  Each click will double the size of the 
	plot.  Double click to return to full scale.

	See also COH, ASD, CRSD, and TFPLOT.

	1994 by Joseph C. Slater
	Modifications:
	--------------
	7/6/00: Changed default FRF calculation from H2 to H1
		 Added H1, H2, and Hv options.
	4/13/15: Converted to Python
	"""

	if len(dt)>1:
		dt=dt(2)-dt(1)

	#if frftype==1
	## freq,Pff=asd(f,dt);
	## freq,Pxf=crsd(x,f,dt);
	## tfunc=Pxf/Pff;
	## elseif frftype==2
	freq,Pxx=asd(x,dt)
	freq,Pxf=crsd(x,f,dt)
	tfunc=Pxx/Pxf
	## elseif frftype==3
	##   [freq,Pxx]=asd(x,dt,n,'noave');
	##   [freq,Pff]=asd(f,dt,n,'noave');
	##   [freq,Pxf]=crsd(x,f,dt,n,'noave');
	##   for i=1:size(Pxx,1)
	##     for j=1:size(Pxx,2)
	##       frfm=[Pff(i,j) Pxf(i,j);Pxf(i,j) Pxx(i,j)];
	##       [v,d]=eig(frfm);
	##       [y,yi]=sort(diag(d));
	##       tfunc(i,j)=-v(yi(1),1)/v(yi(1),2);		
	##     end
	##   end
	## end

	tfunc=np.conj(tfunc)

	## sfft=size(Pxf);
	## if sfft(2)~=1 & ~strcmp(ave,'no')
	##   tfunc=mean(tfunc')';
	## end

	## % If no left hand arguments then plot results
	## if nargout==0
	##   subplot(211)
	##   plot(freq,20*log10(abs(tfunc)))

	##   xlabel('Frequency (Hz)')
	##   ylabel('Mag (dB)')
	##   grid
	##   zoom on
	##   subplot(212)
	##   phase=unwrap(angle(tfunc))*180/pi;
	##   %phase=angle(tfunc)*180/pi;
	##   plot(freq,phase)

	##   xlabel('Frequency (Hz)')
	##   ylabel('Phase (deg)')
	##   grid
	##   phmin_max=[floor(min(phase)/45)*45 ceil(max(phase)/45)*45];
	##   set(gca,'YLim',phmin_max)
	##   gridmin_max=round(phmin_max/90)*90;
	##   set(gca,'YTick',gridmin_max(1):90:gridmin_max(2))
	##   % set(gca,'GridLineStyle','--')
	##   % gridmin_max=round(phmin_max/45)*45;
	##   % set(gca,'YTick',gridmin_max(1):45:gridmin_max(2))
	##   set(gca,'GridLineStyle',':')
	##   set(gca,'YTickLabel',gridmin_max(1):90:gridmin_max(2))
	##   zoom on
	##   return
	## end

	## freqout=freq;
	## tfout=tfunc;
	return freq,tfunc

def frfestHv(x,f,dt):#,n,options)
	"""
	f,Txf = tfestHv(x,f,dt) estimates the 
	Hv Frequency Response Function (FRF) between x and f.
	tfestHv calculates the Hv FRF (Sxf/Sff) by default.
	Freq is the frequency vector in Hertz and Txf is the
	Transfer Function in complex form.
	dt is the time step of the sampled data. 
	If dt is replaced by the time vector, dt will be extracted 
	from the time vector using dt = t[2]-t[1].

	NOT ENABLED YET
	If X and 
	F are matrices, TFEST  will find the frequency response
	function for each column and average the results unless 
	OPTIONS{1} is set to 'no'. 

	NOT ENABLED YET
	TFEST(X,Y,DT,N,OPTIONS) plots the Frequency Response Function 
	if there are no output arguments.  Click in the region of 
	interest to zoom in.  Each click will double the size of the 
	plot.  Double click to return to full scale.

	See also COH, ASD, CRSD, and TFPLOT.

	1994 by Joseph C. Slater
	Modifications:
	--------------
	7/6/00: Changed default FRF calculation from H2 to H1
		 Added H1, H2, and Hv options.
	4/13/15: Converted to Python
	"""

	if len(dt)>1:
		dt=dt(2)-dt(1)

	#if frftype==1
	## freq,Pff=asd(f,dt);
	## freq,Pxf=crsd(x,f,dt);
	## tfunc=Pxf/Pff;
	## elseif frftype==2
	##   [freq,Pxx]=asd(x,dt,n,'noave');
	##   [freq,Pxf]=crsd(x,f,dt,n,'noave');
	##   tfunc=Pxx./Pxf;
	## elseif frftype==3
	freq,Pxx=asd(x,dt)
	freq,Pff=asd(f,dt)
	freq,Pxf=crsd(x,f,dt)
	for i in arange(0,len(Pxx)):
	##	for j in arange(0,len(Pxx)):
		frfm=array([Pff[i,0], Pxf[i,0]],[Pxf[i,0], Pxx[i,0]])
		lam,vecs=la.eigh(frfm)
		index=lam.argsort()
		lam=lam[index]
		vecs=vecs[:,index]
		tfunc[i]=-vec[0,0]/v[0,1]		


	tfunc=np.conj(tfunc)

	## sfft=size(Pxf);
	## if sfft(2)~=1 & ~strcmp(ave,'no')
	##   tfunc=mean(tfunc')';
	## end

	## % If no left hand arguments then plot results
	## if nargout==0
	##   subplot(211)
	##   plot(freq,20*log10(abs(tfunc)))

	##   xlabel('Frequency (Hz)')
	##   ylabel('Mag (dB)')
	##   grid
	##   zoom on
	##   subplot(212)
	##   phase=unwrap(angle(tfunc))*180/pi;
	##   %phase=angle(tfunc)*180/pi;
	##   plot(freq,phase)

	##   xlabel('Frequency (Hz)')
	##   ylabel('Phase (deg)')
	##   grid
	##   phmin_max=[floor(min(phase)/45)*45 ceil(max(phase)/45)*45];
	##   set(gca,'YLim',phmin_max)
	##   gridmin_max=round(phmin_max/90)*90;
	##   set(gca,'YTick',gridmin_max(1):90:gridmin_max(2))
	##   % set(gca,'GridLineStyle','--')
	##   % gridmin_max=round(phmin_max/45)*45;
	##   % set(gca,'YTick',gridmin_max(1):45:gridmin_max(2))
	##   set(gca,'GridLineStyle',':')
	##   set(gca,'YTickLabel',gridmin_max(1):90:gridmin_max(2))
	##   zoom on
	##   return
	## end

	## freqout=freq;
	## tfout=tfunc;
	return freq,tfunc






	## def acorr(self, x, **kwargs):
    ##     """
    ##     Plot the autocorrelation of `x`.

    ##     Parameters
    ##     ----------

    ##     x : sequence of scalar

    ##     hold : boolean, optional, default: True

    ##     detrend : callable, optional, default: `mlab.detrend_none`
    ##         x is detrended by the `detrend` callable. Default is no
    ##         normalization.

    ##     normed : boolean, optional, default: True
    ##         if True, normalize the data by the autocorrelation at the 0-th
    ##         lag.

    ##     usevlines : boolean, optional, default: True
    ##         if True, Axes.vlines is used to plot the vertical lines from the
    ##         origin to the acorr. Otherwise, Axes.plot is used.

    ##     maxlags : integer, optional, default: 10
    ##         number of lags to show. If None, will return all 2 * len(x) - 1
    ##         lags.

    ##     Returns
    ##     -------
    ##     (lags, c, line, b) : where:

    ##       - `lags` are a length 2`maxlags+1 lag vector.
    ##       - `c` is the 2`maxlags+1 auto correlation vectorI
    ##       - `line` is a `~matplotlib.lines.Line2D` instance returned by
    ##         `plot`.
    ##       - `b` is the x-axis.

    ##     Other parameters
    ##     -----------------
    ##     linestyle : `~matplotlib.lines.Line2D` prop, optional, default: None
    ##         Only used if usevlines is False.

    ##     marker : string, optional, default: 'o'

    ##     Notes
    ##     -----
    ##     The cross correlation is performed with :func:`numpy.correlate` with
    ##     `mode` = 2.

    ##     Examples
    ##     --------

    ##     `~matplotlib.pyplot.xcorr` is top graph, and
    ##     `~matplotlib.pyplot.acorr` is bottom graph.

    ##     .. plot:: mpl_examples/pylab_examples/xcorr_demo.py

    ##     """
    ##     return self.xcorr(x, x, **kwargs)

    ## @docstring.dedent_interpd
    ## def xcorr(self, x, y, normed=True, detrend=mlab.detrend_none,
    ##           usevlines=True, maxlags=10, **kwargs):
    ##     """
    ##     Plot the cross correlation between *x* and *y*.

    ##     Parameters
    ##     ----------

    ##     x : sequence of scalars of length n

    ##     y : sequence of scalars of length n

    ##     hold : boolean, optional, default: True

    ##     detrend : callable, optional, default: `mlab.detrend_none`
    ##         x is detrended by the `detrend` callable. Default is no
    ##         normalization.

    ##     normed : boolean, optional, default: True
    ##         if True, normalize the data by the autocorrelation at the 0-th
    ##         lag.

    ##     usevlines : boolean, optional, default: True
    ##         if True, Axes.vlines is used to plot the vertical lines from the
    ##         origin to the acorr. Otherwise, Axes.plot is used.

    ##     maxlags : integer, optional, default: 10
    ##         number of lags to show. If None, will return all 2 * len(x) - 1
    ##         lags.

    ##     Returns
    ##     -------
    ##     (lags, c, line, b) : where:

    ##       - `lags` are a length 2`maxlags+1 lag vector.
    ##       - `c` is the 2`maxlags+1 auto correlation vectorI
    ##       - `line` is a `~matplotlib.lines.Line2D` instance returned by
    ##         `plot`.
    ##       - `b` is the x-axis (none, if plot is used).

    ##     Other parameters
    ##     -----------------
    ##     linestyle : `~matplotlib.lines.Line2D` prop, optional, default: None
    ##         Only used if usevlines is False.

    ##     marker : string, optional, default: 'o'

    ##     Notes
    ##     -----
    ##     The cross correlation is performed with :func:`numpy.correlate` with
    ##     `mode` = 2.
    ##     """

    ##     Nx = len(x)
    ##     if Nx != len(y):
    ##         raise ValueError('x and y must be equal length')

    ##     x = detrend(np.asarray(x))
    ##     y = detrend(np.asarray(y))

    ##     c = np.correlate(x, y, mode=2)

    ##     if normed:
    ##         c /= np.sqrt(np.dot(x, x) * np.dot(y, y))

    ##     if maxlags is None:
    ##         maxlags = Nx - 1

    ##     if maxlags >= Nx or maxlags < 1:
    ##         raise ValueError('maglags must be None or strictly '
    ##                          'positive < %d' % Nx)

    ##     lags = np.arange(-maxlags, maxlags + 1)
    ##     c = c[Nx - 1 - maxlags:Nx + maxlags]

    ##     if usevlines:
    ##         a = self.vlines(lags, [0], c, **kwargs)
    ##         b = self.axhline(**kwargs)
    ##     else:

    ##         kwargs.setdefault('marker', 'o')
    ##         kwargs.setdefault('linestyle', 'None')
    ##         a, = self.plot(lags, c, **kwargs)
    ##         b = None
    ##     return lags, c, a, b
if __name__ == "__main__":
    import doctest
    doctest.testmod()