Python TimeSeries.peek方法代码示例

# 需要导入模块: from timeseries import TimeSeries [as 别名]
# 或者: from timeseries.TimeSeries import peek [as 别名]

#.........这里部分代码省略.........
                # Set the scale type on each axis
                ax.set_xscale('log')

                # Set the formatting of the tick labels
                xformatter = plt.FuncFormatter(log_10_product)
                ax.xaxis.set_major_formatter(xformatter)

                # Geometric mean
                ax.plot(freqs, logiobs / np.log(10.0),  color='k', label='geometric mean of power spectra at each pixel')
                #ax.plot(freqs, bf2, color='k', label='best fit n=%4.2f +/- %4.2f' % (param2[1], nerr2))

                # Power at each frequency - distributions
                ax.plot(freqs, np.log10(lim[0, 0, :]), label=s_L68.label, color=s_L68.color, linewidth=s_L68.linewidth, linestyle=s_L68.linestyle)
                ax.plot(freqs, np.log10(lim[1, 0, :]), label=s_L95.label, color=s_L95.color, linewidth=s_L95.linewidth, linestyle=s_L95.linestyle)
                ax.plot(freqs, np.log10(lim[0, 1, :]), label=s_U68.label, color=s_U68.color, linewidth=s_U68.linewidth, linestyle=s_U68.linestyle)
                ax.plot(freqs, np.log10(lim[1, 1, :]), label=s_U95.label, color=s_U95.color, linewidth=s_U95.linewidth, linestyle=s_U95.linestyle)

                # Position of the fitted peak in each distribution
                ax.plot(freqs, logiobs_peak_location / np.log(10.0),  color='m', label='fitted frequency')

                # Extra information for the plot
                ax.axvline(five_min, color=s5min.color, linestyle=s5min.linestyle, label=s5min.label)
                ax.axvline(three_min, color=s3min.color, linestyle=s3min.linestyle, label=s3min.label)
                plt.xlabel('frequency (%s)' % (freqfactor[1]))
                plt.ylabel('power [%i time series, %i samples each]' % (nx * ny, nt))
                plt.title(data_name + ' : geometric mean')
                plt.legend(loc=3, fontsize=10, framealpha=0.5)
                plt.savefig(savefig + '.geometric_mean_power_spectra.%s' % (savefig_format))
                plt.close('all')
                # -------------------------------------------------------------

                # plot out the time series
                plt.figure(4)
                full_ts.peek()
                plt.savefig(savefig + '.full_ts_timeseries.%s' % (savefig_format))
                plt.close('all')

                # -------------------------------------------------------------
                # plot some histograms of the power at a small number of
                # frequencies.
                """
                histogram_loc2, hpwr2, lim2 = calculate_histograms(nposfreq, pwr, 100)

                findex = []
                f_of_interest = [0.5 * five_min, five_min, three_min, 2 * three_min, 3 * three_min]
                for thisf in f_of_interest:
                    findex.append(np.unravel_index(np.argmin(np.abs(thisf - freqs)), freqs.shape)[0])
                plt.figure(3)
                plt.xlabel('power')
                plt.ylabel('proportion found at given frequency')
                plt.title(data_name + ' - power distributions')
                for f in findex:
                    xx = histogram_loc2[1:] / np.log(10.0)
                    yy = hpwr2[f, :]
                    plt.loglog(xx, yy, label='%7.2f %s' % (freqs[f], freqfactor[1]))
                plt.legend(loc=3, fontsize=10, framealpha=0.5)
                plt.savefig(savefig + '.notlog_power_spectra_distributions.%s' % (savefig_format))

                # plot out the time series
                plt.figure(4)
                full_ts.peek()
                plt.savefig(savefig + '.full_ts_timeseries.%s' % (savefig_format))
                plt.close('all')
                """
                ###############################################################
                # Time series plots

# 需要导入模块: from timeseries import TimeSeries [as 别名]
# 或者: from timeseries.TimeSeries import peek [as 别名]

#.........这里部分代码省略.........
                        lim[i, 1, f] = np.exp(h[1][hi])

                # Give the best plot we can under the circumstances.  Since we have been
                # looking at the log of the power, plots are slightly different
                plt.figure(2)
                plt.loglog(freqs, np.exp(logiobs), label='geometric mean of power spectra at each pixel')
                plt.loglog(freqs, bf2, color='k', label='best fit n=%4.2f +/- %4.2f' % (param2[1], nerr2))
                plt.loglog(freqs, lim[0, 0, :], linestyle='--', label='lower 68%')
                plt.loglog(freqs, lim[0, 1, :], linestyle='--', label='upper 68%')
                plt.loglog(freqs, lim[1, 0, :], linestyle=':', label='lower 95%')
                plt.loglog(freqs, lim[1, 1, :], linestyle=':', label='upper 95%')
                plt.axvline(five_min, color='k', linestyle='-.', label='5 mins.')
                plt.axvline(three_min, color='k', linestyle='--', label='3 mins.')
                plt.xlabel('frequency (Hz)')
                plt.ylabel('power [%i time series, %i samples each]' % (nx * ny, nt))
                plt.title(data_name + ' - gPS')
                plt.legend(loc=1, fontsize=10)
                plt.savefig(savefig + '.geometric_mean_power_spectra.png')

                # plot some histograms of the log power at a small number of equally spaced
                # frequencies
                findex = [0, 11, 19, 38, 76]
                plt.figure(3)
                plt.xlabel('$\log_{10}(power)$')
                plt.ylabel('proportion found at given frequency')
                plt.title(data_name + ' - power distributions')
                for f in findex:
                    plt.plot(h[1][1:] / np.log(10.0), hpwr[f, :], label='%7.5f Hz' % (freqs[f]))
                plt.legend(loc=3, fontsize=10)
                plt.savefig(savefig + '.power_spectra_distributions.png')

                # plot out the time series
                plt.figure(4)
                full_ts.peek()
                plt.savefig(savefig + '.full_ts_timeseries.png')
                plt.close('all')

                ###############################################################
                # Time series plots
                # Plot all the analyzed time series
                plt.figure(10)
                for i in range(0, nx):
                    for j in range(0, ny):
                        plt.plot(t, dc_analysed[j, i, :])
                plt.xlabel('time (seconds)')
                plt.ylabel('analyzed emission ' + tsdetails)
                plt.title(data_name)
                plt.ylim(dc_analysed_minmax)
                plt.xlim((t[0], t[-1]))
                plt.savefig(savefig + '.all_analyzed_ts.png')

                # Plot a histogram of the studied data at each time
                bins = 50
                hist_dc_analysed = np.zeros((bins, nt))
                for this_time in range(0, nt):
                    hist_dc_analysed[:, this_time], bin_edges = np.histogram(dc_analysed[:, :, this_time], bins=bins, range=dc_analysed_minmax)
                hist_dc_analysed = hist_dc_analysed / (1.0 * nx * ny)
                plt.figure(12)
                plt.xlabel('time (seconds)')
                plt.ylabel('analyzed emission ' + tsdetails)
                plt.imshow(hist_dc_analysed, aspect='auto', origin='lower',
                           extent=(t[0], t[-1], dc_analysed_minmax[0], dc_analysed_minmax[1]))
                plt.colorbar()
                plt.title(data_name)
                plt.savefig(savefig + '.all_analyzed_ts_histogram.png')

# 需要导入模块: from timeseries import TimeSeries [as 别名]
# 或者: from timeseries.TimeSeries import peek [as 别名]
alpha = 0.0001

data[0] = 1.0
for i in range(0, nt - 1):
    data[i+1] = data[i] + alpha*np.random.normal()


ts = TimeSeries(dt * np.arange(0, nt), data)

plt.figure(1)
ts.peek_ps()
plt.loglog()


plt.figure(2)
ts.peek()

this = ([ts.pfreq, ts.ppower],)

norm_estimate = np.zeros((3,))
norm_estimate[0] = ts.ppower[0]
norm_estimate[1] = norm_estimate[0] / 1000.0
norm_estimate[2] = norm_estimate[0] * 1000.0

background_estimate = np.zeros_like(norm_estimate)
background_estimate[0] = np.mean(ts.ppower[-10:-1])
background_estimate[1] = background_estimate[0] / 1000.0
background_estimate[2] = background_estimate[0] * 1000.0

estimate = {"norm_estimate": norm_estimate,
            "background_estimate": background_estimate}

# 需要导入模块: from timeseries import TimeSeries [as 别名]
# 或者: from timeseries.TimeSeries import peek [as 别名]
# plot some histograms of the log power at a small number of equally spaced
# frequencies
findex = np.arange(0, nposfreq, nposfreq / 5)
plt.figure(3)
plt.xlabel('$\log_{10}(power)$')
plt.ylabel('proportion found at given frequency')
plt.title(data_name + ' - power distributions')
for f in findex:
    plt.plot(h[1][1:] / np.log(10.0), hpwr[f, :], label='%7.5f Hz' % (freqs[f]))
plt.legend(loc=3, fontsize=10)
plt.savefig(savefig + '.power_spectra_distributions.png')


# plot out the time series
plt.figure(4)
full_ts.peek()
plt.savefig(savefig + '.full_ts_timeseries.png')



#
# Make maps of the Fourier power
#
fmap = []
franges = [[1.0/360.0, 1.0/240.0], [1.0/240.0, 1.0/120.0]]
for fr in franges:
    ind = []
    for i, testf in enumerate(freqs):
        if testf >= fr[0] and testf <= fr[1]:
            ind.append(i)
    fmap.append(np.sum(pwr[:,:,ind[:]], axis=2))

# 需要导入模块: from timeseries import TimeSeries [as 别名]
# 或者: from timeseries.TimeSeries import peek [as 别名]
            raise ValueError, "Input vector needs to be bigger than window size."
    if window_len < 3:
            return x
    if not window in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']:
            raise ValueError, "Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'"
    s = np.r_[2 * x[0] - x[window_len - 1::-1], x, 2 * x[-1] - x[-1:-window_len:-1]]
    if window == 'flat': #moving average
            w = np.ones(window_len, 'd')
    else:
            w = eval('np.' + window + '(window_len)')
    y = np.convolve(w / w.sum(), s, mode='same')
    return y[window_len:-window_len + 1]

tsoriginal = TimeSeries(t, data)
plt.figure(10)
tsoriginal.peek()

meandata = np.mean(data)

# relative
data = (data - meandata) / meandata


#data = data - smooth(data, window_len=84)

# Create a time series object
ts = TimeSeries(t, data)
ts.label = 'emission'
ts.units = 'arb. units'
ts.name = 'simulated data [n=%4.2f]' % (model_param[1])