def plot_psd_resampled():
    global t, ch1, ch2, sps, sampl, title
    
    plt.clf()
    nfft=2048 #int(2**np.ceil(np.log2(sampl)));

    ax1=plt.subplot(211)
    (Pxx, freqs)=plt.psd(ch1[::5], NFFT=nfft, Fs=sps/5, window=mlab.window_hanning, color='red')
    plt.title('PSD (NFFT='+str(nfft)+', Fs='+str(sps/5)+'Hz, hann) from '+title)
    plt.xlabel('')

    ax2=plt.subplot(212)
    (Pxx, freqs)=plt.psd(ch2[::5], NFFT=nfft, Fs=sps/5, window=mlab.window_hanning, color='green')
    plt.ylabel('')
    ax2.set_ylim(ax1.get_ylim())

    fig=plt.gcf();
    f_dpi=fig.get_dpi()
    #print "dpi: %i" % (f_dpi)
    
    #f_size=fig.get_size_inches()
    #print "figure size in Inches", f_size
    #print "or pixels %i x %i pixel" % (f_dpi*f_size[0], f_dpi*f_size[1])
    #fig.set_size_inches([f_size[0]*2, f_size[1]]);
    
    fig.set_size_inches(1200.0/f_dpi, 500.0/f_dpi);
    
    #f_size=fig.get_size_inches()
    #print "resized to %i x %i pixel" % (f_dpi*f_size[0], f_dpi*f_size[1])

    
    plt.savefig(os.path.join(_LOCATION_CHARTS,"last_psd_res.png"), dpi=96, bbox_inches='tight')
    plt.clf()
    print "plotting psd_resampled finished"

def fftplt1(x,Mch,dtrnd=True,demean=True):
	fs=1

	if demean:
		x = x - x.mean()

	if dtrnd:
		Pout = plt.psd(x, NFFT=Mch, detrend=pylab.detrend_linear, noverlap=Mch/2, Fs=fs)
	else:
		Pout = plt.psd(x, NFFT=Mch, detrend=pylab.detrend_none, noverlap=Mch/2, Fs=fs)

	xdtrnd = pylab.detrend_linear(x)
	xauto = mcmath.acorr_mlab(xdtrnd,2)
	rhoxauto = (xauto[1] + np.sqrt(abs(xauto[2])))/2
	R = mcmath.redspec(rhoxauto,np.arange(Mch/2),Mch)

	P = Pout[0][:R.size]
	F = Pout[1][:R.size]
	Psum = P.sum()
	Rsum = R.sum()
	PRratio = Psum/Rsum
	Rcmp = R*PRratio

	plt.figure()
	plt.plot(F,P)
	plt.plot(F,Rcmp)

	return (F,P,Rcmp)

 def inspect_source_psd(self, ic):
     data = self._get_ica_data()
     source = self.ica._transform_epochs(data, concatenate=True)[ic]
     sfreq = data.info['sfreq']
     plt.figure()
     plt.psd(source, Fs=sfreq, NFFT=128, noverlap=0, pad_to=None)
     plt.show()

def fftplt3(x,Mch,pval=0.1,dtrnd=True,demean=True,titlestr='',Dt=0.01):
	"""Differs from fftplt3 solely by plotting the power spectral density
	in log-form, 10*log10(P)."""
	fs=1
	titl_plc = (1.,1.)

	if demean:
		x = x - x.mean()

	if dtrnd:
		Pout = plt.psd(x, NFFT=Mch, detrend=pylab.detrend_linear, noverlap=Mch/2, Fs=fs)
	else:
		Pout = plt.psd(x, NFFT=Mch, detrend=pylab.detrend_none, noverlap=Mch/2, Fs=fs)

	xdtrnd = pylab.detrend_linear(x)
	xauto = mcmath.acorr_mlab(xdtrnd,2)
	rhoxauto = (xauto[1] + np.sqrt(abs(xauto[2])))/2
	R = mcmath.redspec(rhoxauto,np.arange(Mch/2),Mch)

	P = Pout[0][:R.size]
	F = Pout[1][:R.size]
	Psum = P.sum()
	Rsum = R.sum()
	PRratio = Psum/Rsum
	Rcmp = R*PRratio

	dof = (x.size / (Mch/2.)) * 1.2
	Fval = stats.f.isf(pval,dof,dof)

	tst = P / Rcmp
	pass1 = np.where(tst > Fval)[0]
	maxs = mcmath.extrema_find(P,'max',t=F,dt=Dt)
	max_ind = np.intersect1d(maxs,pass1)
	Fmaxs = F[max_ind]
	Fmaxs2 = np.append(Fmaxs,0.1)
	Fmaxs2.sort()
	Tmaxs = 1/(Fmaxs2)
	Tmaxs = np.round(Tmaxs,2)

	ax = plt.gca()
	ax.plot(F,10*np.log10(Rcmp))

	ax.set_xticks(Fmaxs2)
	ax.set_xticklabels(Tmaxs)
	my_fmt_xdate(ax,rot=90,hal='center')
	multiline(Fmaxs,c='red',ls='--')

	xtl = ax.get_xticklabels()
	ytl = ax.get_yticklabels()
	plt.setp(xtl,'size',10)
	plt.setp(ytl,'size',9)

	ppct = int((1 - pval)*100)
	titl_str = '%s FFT (chunksize: %d) :: peak CL: %d%%' % (titlestr,Mch,ppct)
	plt.text(titl_plc[0],titl_plc[1],titl_str,ha='right',va='bottom',size=12,transform=ax.transAxes)
	plt.xlabel('Peak period (yrs)',size=11)


	return (F,P,Rcmp)

def do_check_spectrum(hostData, DUTData, samplingRate, fLow, fHigh, margainLow, margainHigh):
    # reduce FFT resolution to have averaging effects
    N = 512 if (len(hostData) > 512) else len(hostData)
    iLow = N * fLow / samplingRate + 1  # 1 for DC
    if iLow > (N / 2 - 1):
        iLow = N / 2 - 1
    iHigh = N * fHigh / samplingRate + 1  # 1 for DC
    if iHigh > (N / 2 + 1):
        iHigh = N / 2 + 1
    print fLow, iLow, fHigh, iHigh, samplingRate

    Phh, freqs = plt.psd(
        hostData,
        NFFT=N,
        Fs=samplingRate,
        Fc=0,
        detrend=plt.mlab.detrend_none,
        window=plt.mlab.window_hanning,
        noverlap=0,
        pad_to=None,
        sides="onesided",
        scale_by_freq=False,
    )
    Pdd, freqs = plt.psd(
        DUTData,
        NFFT=N,
        Fs=samplingRate,
        Fc=0,
        detrend=plt.mlab.detrend_none,
        window=plt.mlab.window_hanning,
        noverlap=0,
        pad_to=None,
        sides="onesided",
        scale_by_freq=False,
    )
    print len(Phh), len(Pdd)
    print "Phh", abs(Phh[iLow:iHigh])
    print "Pdd", abs(Pdd[iLow:iHigh])
    amplitudeRatio = np.sqrt(abs(Pdd[iLow:iHigh] / Phh[iLow:iHigh]))
    ratioMean = np.mean(amplitudeRatio)
    amplitudeRatio = amplitudeRatio / ratioMean
    print "Normialized ratio", amplitudeRatio
    print "ratio mean for normalization", ratioMean
    positiveMax = abs(max(amplitudeRatio))
    negativeMin = abs(min(amplitudeRatio))
    passFail = (
        True if (positiveMax < (margainHigh / 100.0 + 1.0)) and ((1.0 - negativeMin) < margainLow / 100.0) else False
    )
    RatioResult = np.zeros(len(amplitudeRatio), dtype=np.int16)
    for i in range(len(amplitudeRatio)):
        RatioResult[i] = amplitudeRatio[i] * 1024  # make fixed point
    print "positiveMax", positiveMax, "negativeMin", negativeMin
    return (passFail, negativeMin, positiveMax, RatioResult)

def problem_5():
    """
    This function uses pyplot.psd to plot the spectrum for each of the time
    series generated above.

    Inputs:
        None

    Outputs:
        This will automatically generate a plot for each of the 4 data-sets.
    """

    plt.figure()
    plt.psd(gdp)
    plt.title('GDP Spectrum')
    plt.show()

    plt.figure()
    plt.psd(cpi)
    plt.title('CPI Spectrum')
    plt.show()

    plt.figure()
    plt.psd(cons)
    plt.title('Consumption Spectrum')
    plt.show()

    plt.figure()
    plt.psd(inv)
    plt.title('Investment Spectrum')
    plt.show()

def chickling_csd_2ch(shotno, date=time.strftime("%Y%m%d")):

	fname, data = file_finder(shotno,date)
	data = quickextract_data(fname)

	#reshape the array of x points (20M for 1s) into a 2d array each with 40k segments.
	phasediff_co2 = np.unwrap(data[0]['phasediff_co2'])
	phasediff_hene = np.unwrap(data[0]['phasediff_hene'])
	fs = data[0]['samplerate']
	
	plt.figure("2 Channels | Blue = PSD Scene | Orange = PSD Reference | Green = CSD | shot " + str(shotno) +  " Date " + str(date))
	plt.psd(phasediff_co2, Fs=fs)
	plt.psd(phasediff_hene,Fs=fs)
	plt.csd(phasediff_co2, phasediff_hene, Fs=fs)
	plt.show()

def compute_and_plot_psd(data, sfreq, NFFT=512, show=True):
    """ Computes the power spectral density and produces a plot.

    Parameters
    ----------
    data : ndarray (n_times)
        The signal.
    sfreq : float
        Sampling frequency.
    NFFT : int (power of 2)
        Number of bins for each block of FFT.
    show : bool
        Display or hide plot. (Default is True)

    Returns
    -------
    power : ndarray
        The power spectral density of the signal.
    freqs : ndarray
        Frequencies.
    """

    if show is False:
        pl.ioff()
    power, freqs = pl.psd(data, Fs=sfreq, NFFT=NFFT)

    return power, freqs

def calculate_amplitude(st):
    """
    Calculate a noise proxy for the first trace in that station.
    """
    res = plt.psd(st[0].data,Fs=st[0].stats.sampling_rate,NFFT=4096,noverlap=4096/2)
    amp_sum = np.sum(res[0])
    return amp_sum

def _plot_fft(tdms, **kargs):
	ti = kargs['ti']
	fft_len = kargs['fft_len']
	fft_scale = kargs['fft_scale']
	p, f = plt.psd(tdms.wav.__getslice__(*ti),
			fft_len,
			tdms.fs)
	if fft_scale == 'db':
		plt.ylabel(r'Amplitude $\frac{mV^2}{Hz}$ ($dB$)')
	elif fft_scale == 'linear':
		plt.cla()
		plt.plot(f, p)
		plt.ylabel(r'Amplitude $\frac{mV^2}{Hz}$')
	plt.xlabel(u'Frequency ($Hz$)')
	return
	fft = np.fft.fft(tdms.wav.__getslice__(*ti), fft_len)[:fft_len/2]

	if fft_unit == 'pow':
		fft = fft * fft.conjugate() / fft_len
	elif fft_unit == 'amp':
		fft = np.abs(fft / fft_len)
		plt.ylabel('Amplitude mV/Hz')

	plt.plot(np.fft.fftfreq(fft_len, 1/tdms.fs)[:fft_len/2], fft)
	plt.grid()

def do_check_spectrum_playback(hostData, samplingRate, fLow, fHigh, margainLow, margainHigh):
    # reduce FFT resolution to have averaging effects
    N = 512 if (len(hostData) > 512) else len(hostData)
    iLow = N * fLow / samplingRate + 1 # 1 for DC
    if iLow > (N / 2 - 1):
        iLow = (N / 2 - 1)
    iHigh = N * fHigh / samplingRate + 1 # 1 for DC
    if iHigh > (N / 2 + 1):
        iHigh = N / 2 + 1
    print fLow, iLow, fHigh, iHigh, samplingRate

    Phh, freqs = plt.psd(hostData, NFFT=N, Fs=samplingRate, Fc=0, detrend=plt.mlab.detrend_none,\
        window=plt.mlab.window_hanning, noverlap=0, pad_to=None, sides='onesided',\
        scale_by_freq=False)
    print len(Phh)
    print "Phh", abs(Phh[iLow:iHigh])
    spectrum = np.sqrt(abs(Phh[iLow:iHigh]))
    spectrumMean = np.mean(spectrum)
    spectrum = spectrum / spectrumMean
    print "Mean ", spectrumMean
    print "Normalized spectrum", spectrum
    positiveMax = abs(max(spectrum))
    negativeMin = abs(min(spectrum))
    passFail = True if (positiveMax < (margainHigh / 100.0 + 1.0)) and\
        ((1.0 - negativeMin) < margainLow / 100.0) else False
    spectrumResult = np.zeros(len(spectrum), dtype=np.int16)
    for i in range(len(spectrum)):
        spectrumResult[i] = spectrum[i] * 1024 # make fixed point
    print "positiveMax", positiveMax, "negativeMin", negativeMin
    return (passFail, negativeMin, positiveMax, spectrumResult)

def plot_psd():
    global t, ch1, ch2, sps, sampl, title
    
    plt.clf()
    nfft=2048 #int(2**np.ceil(np.log2(sampl)));

    ax1=plt.subplot(211)
    (Pxx, freqs)=plt.psd(ch1, NFFT=nfft, Fs=sps, window=mlab.window_hanning)
    plt.title('PSD (NFFT='+str(nfft)+', Fs='+str(sps)+'Hz, hann) from '+title)
    plt.xlabel('')

    ax2=plt.subplot(212)
    (Pxx, freqs)=plt.psd(ch2, NFFT=nfft, Fs=sps, window=mlab.window_hanning)
    plt.ylabel('')
    ax2.set_ylim(ax1.get_ylim())

    fig=plt.gcf()
    f_dpi=fig.get_dpi()    
    fig.set_size_inches(1200.0/f_dpi, 500.0/f_dpi);

    plt.savefig(os.path.join(_LOCATION_CHARTS,"last_psd.png"), dpi=96, bbox_inches='tight')
    plt.clf()
    print "plotting psd finished"

def plot_psd_OFDM_symbols():  # pragma: no cover
    """Plot the power spectral density of OFDM modulated symbols.

    This function is not used in any unittest, but it is interesting to
    visualize that the modulate method of the OFDM class is working as it
    should.
    """
    from matplotlib import pyplot as plt
    # xxxxxxxxxx OFDM Details xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
    fft_size = 64
    cp_size = 12
    num_used_subcarriers = 52
    ofdm_object = ofdm.OFDM(fft_size, cp_size, num_used_subcarriers)
    # xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

    # xxxxxxxxxx Input generation (not part of OFDM) xxxxxxxxxxxxxxxxxxxxxx
    num_bits = 2500
    # generating 1's and 0's
    ip_bits = np.random.random_integers(0, 1, num_bits)
    # Number of modulated symbols
    # num_mod_symbols = num_bits * 1
    # xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

    # BPSK modulation
    # bit0 --> -1
    # bit1 --> +1
    ip_mod = 2 * ip_bits - 1

    # OFDM Modulation
    output = ofdm_object.modulate(ip_mod)

    # xxxxxxxxxx Plot the PSD xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
    # MATLAB code to plot the power spectral density
    # close all
    fsMHz = 20e6
    Pxx, W = plt.psd(output, NFFT=fft_size, Fs=fsMHz)
    # [Pxx,W] = pwelch(output,[],[],4096,20);
    plt.plot(
        W,
        10 * np.log10(Pxx)
    )
    plt.xlabel('frequency, MHz')
    plt.ylabel('power spectral density')
    plt.title('Transmit spectrum OFDM (based on 802.11a)')
    plt.show()

def test_compares_psd():
    """Test PSD estimation on raw for plt.psd and scipy.signal.welch
    """
    raw = io.read_raw_fif(raw_fname)

    exclude = raw.info['bads'] + ['MEG 2443', 'EEG 053']  # bads + 2 more

    # picks MEG gradiometers
    picks = pick_types(raw.info, meg='grad', eeg=False, stim=False,
                       exclude=exclude)[:2]

    tmin, tmax = 0, 10  # use the first 60s of data
    fmin, fmax = 2, 70  # look at frequencies between 5 and 70Hz
    n_fft = 2048

    # Compute psds with the new implementation using Welch
    psds_welch, freqs_welch = psd_welch(raw, tmin=tmin, tmax=tmax,
                                        fmin=fmin, fmax=fmax,
                                        proj=False, picks=picks,
                                        n_fft=n_fft, n_jobs=1)

    # Compute psds with plt.psd
    start, stop = raw.time_as_index([tmin, tmax])
    data, times = raw[picks, start:(stop + 1)]
    from matplotlib.pyplot import psd
    out = [psd(d, Fs=raw.info['sfreq'], NFFT=n_fft) for d in data]
    freqs_mpl = out[0][1]
    psds_mpl = np.array([o[0] for o in out])

    mask = (freqs_mpl >= fmin) & (freqs_mpl <= fmax)
    freqs_mpl = freqs_mpl[mask]
    psds_mpl = psds_mpl[:, mask]

    assert_array_almost_equal(psds_welch, psds_mpl)
    assert_array_almost_equal(freqs_welch, freqs_mpl)

    assert_true(psds_welch.shape == (len(picks), len(freqs_welch)))
    assert_true(psds_mpl.shape == (len(picks), len(freqs_mpl)))

    assert_true(np.sum(freqs_welch < 0) == 0)
    assert_true(np.sum(freqs_mpl < 0) == 0)

    assert_true(np.sum(psds_welch < 0) == 0)
    assert_true(np.sum(psds_mpl < 0) == 0)

def _plot_freq_hist(tdms, **kargs):
	ti = kargs['ti']
	fi = kargs['fi']
	fft_len = kargs['fft_len']
	fft_scale = kargs['fft_scale']
	n_bins = kargs['n_bins']

#	bins = list()
	s, f = plt.psd(tdms.wav.__getslice__(*ti),
			fft_len,
			tdms.fs)
	sfi = None
	ffi = None
	for v, i in zip(f, xrange(len(f))):
		if fi[0] is not None and v > fi[0]:
			if sfi is None:
				sfi = i
		if fi[1] is not None and v > fi[1]:
			ffi = i
			break
	if sfi is None:
		sfi = 0
	if ffi is None:
		ffi = len(f)
	s = s[sfi:ffi]
	plt.cla()
#	s = tdms.wav.__getslice__(*ti)
	t = sum(s)
	m = len(s) / n_bins
	r = len(s) - m * n_bins
	for i in xrange(n_bins):
		step = i * m + (i if i < r else r)
		b = step, step + m + (1 if i < r else 0)
#		bins.append(b, sum(s.__getslice__(*b)))
		plt.bar(b[0] + sfi, sum(s.__getslice__(*b))/t, b[1] - b[0])
	plt.ylabel(r'Amplitude $\frac{mV^2}{Hz}$')
	plt.xlabel(u'Frequency ($Hz$)')

raw.info['bads'] = ['MEG 2443', 'EEG 053']  # mark bad channels

# Set up pick list: Gradiometers - bad channels
picks = mne.pick_types(raw.info, meg='grad', exclude='bads')

order = 5  # define model order
picks = picks[:1]

# Estimate AR models on raw data
b, a = fit_iir_model_raw(raw, order=order, picks=picks, tmin=60, tmax=180)
d, times = raw[0, 10000:20000]  # look at one channel from now on
d = d.ravel()  # make flat vector
innovation = signal.convolve(d, a, 'valid')
d_ = signal.lfilter(b, a, innovation)  # regenerate the signal
d_ = np.r_[d_[0] * np.ones(order), d_]  # dummy samples to keep signal length

###############################################################################
# Plot the different time series and PSDs
plt.close('all')
plt.figure()
plt.plot(d[:100], label='signal')
plt.plot(d_[:100], label='regenerated signal')
plt.legend()

plt.figure()
plt.psd(d, Fs=raw.info['sfreq'], NFFT=2048)
plt.psd(innovation, Fs=raw.info['sfreq'], NFFT=2048)
plt.psd(d_, Fs=raw.info['sfreq'], NFFT=2048, linestyle='--')
plt.legend(('Signal', 'Innovation', 'Regenerated signal'))
plt.show()

import argparse
import sys

# Calculates sampling rate and plots 4 channel specgram of EEG data.csv passed in as in command line arg

eegData = pd.read_csv(sys.argv[1])
#eegData = pd.read_csv('Filtered_EEG3.csv')
time = (eegData['Timestamp (ms)'][eegData.shape[0] - 1] - eegData['Timestamp (ms)'][0]) / 1000
samples = eegData.shape[0] - 1
samplingRate = samples/time
print('estimated sampling rate: ',  samplingRate)


plt.figure()
plt.subplot(2,2,1)
plt.psd(eegData['Electrode 1'], Fs=samplingRate)
plt.xlim(0,65)
plt.ylim(-60,20)

plt.subplot(2,2,2)
plt.psd(eegData['Electrode 2'], Fs=samplingRate)
plt.xlim(0,65)
plt.ylim(-60,20)
plt.ylabel('')

plt.subplot(2,2,3)
plt.psd(eegData['Electrode 3'], Fs=samplingRate)
plt.xlim(0,65)
plt.ylim(-60,20)
plt.ylabel('')

def power_spectral_density(data_array, frequency, segment_size=256, window_method=False):
    """Return the power spectral density using matplotlib.pyplot.psd function."""
    return psd(data_array, NFFT=segment_size, Fs = frequency)

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab

fs = 1000
t = np.linspace(0, 0.3, 301)
A = np.array([2, 8]).reshape(-1, 1)
f = np.array([150, 140]).reshape(-1, 1)
xn = (A * np.exp(2j * np.pi * f * t)).sum(axis=0) + 5 * np.random.randn(*t.shape)

yticks = np.arange(-50, 30, 10)
xticks = np.arange(-500, 550, 100)
plt.subplots_adjust(hspace=0.45, wspace=0.3)
ax = plt.subplot(1, 2, 1)

plt.psd(xn, NFFT=301, Fs=fs, window=mlab.window_none, pad_to=1024,
        scale_by_freq=True)
plt.title('Periodogram')
plt.yticks(yticks)
plt.xticks(xticks)
plt.grid(True)
plt.xlim(-500, 500)

plt.subplot(1, 2, 2, sharex=ax, sharey=ax)
plt.psd(xn, NFFT=150, Fs=fs, window=mlab.window_none, noverlap=75, pad_to=512,
        scale_by_freq=True)
plt.title('Welch')
plt.xticks(xticks)
plt.yticks(yticks)
plt.ylabel('')
plt.grid(True)
plt.xlim(-500, 500)

np.random.seed(19680801)


dt = 0.01
t = np.arange(0, 10, dt)
nse = np.random.randn(len(t))
r = np.exp(-t/0.05)

cnse = np.convolve(nse, r)*dt
cnse = cnse[:len(t)]
s = 0.1*np.sin(2*np.pi*t) + cnse

plt.subplot(211)
plt.plot(t, s)
plt.subplot(212)
plt.psd(s, 512, 1/dt)

plt.show()

###############################################################################
# Compare this with the equivalent Matlab code to accomplish the same thing::
#
#     dt = 0.01;
#     t = [0:dt:10];
#     nse = randn(size(t));
#     r = exp(-t/0.05);
#     cnse = conv(nse, r)*dt;
#     cnse = cnse(1:length(t));
#     s = 0.1*sin(2*pi*t) + cnse;
#
#     subplot(211)