Aim: My aim is the calculation and display of a power-law distribution.
Input to Phyton: The input files for Python provided by me are 1-column .txt files. Each .txt file has a column of numeric values extracted from time-series data. My script loads the .txt files into Python to then apply the Welch method on the former. The result is the power-spectrum of the initial time-series data. All 23 subjects (.txt files) are plotted into one image.
From my knowledge, I need to apply the logarithm on both axes of the power-spectrum to calculate the power-law distribution. I have two questions:
1.There seem to be many methods(tools) to calculate the power-law distribution in Python. Based on the fact that I used Welch instead of the FFT, what method would you recommend to now calculate the power-law distribution?
2.Is it possible to directly calculate the power-law distribution based on the time-series data of my input, or do I need to create the code in such a way that it takes the output of welch (the power-spectrum) as input to finally create the power-law distribution? In other words: can I shorten my code and even skip Welch method if my aim is only to display the power-law distribution (and not the power-spectrum)?
Thanks, Philipp
My current code to create the power-spectrum using Welch method is shown below.
# Initialize
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from matplotlib.pyplot import cm
# Numpy.loadtxt – Loads data from a textfile.
# Scipy.signal.welch – Creation of the power-spectrum via welch method. f, Welch creates the ideal frequencies (f, Welch = Power Spectrum or Power Spectral Density)
Subjects = ["Subject1", "Subject2", "Subject3", "Subject4", "Subject5", "Subject7", "Subject8", "Subject9", "Subject10", "Subject11", "Subject12", "Subject13",
"Subject14", "Subject15", "Subject16", "Subject17", "Subject18", "Subject19", "Subject20", "Subject22", "Subject23", "Subject24", "Subject25"]
for Subject in Subjects:
Txt = np.loadtxt("/volumes/SanDisk2/fmri/dataset/processed/extracted_timeseriespython/restingstate/{0}/TimeSeries.SPC.Core_ROI.{0}.txt".format(Subject), comments="#", delimiter=None,
converters=None, skiprows=0, usecols=0, unpack=False, ndmin=0, encoding=None, max_rows=None, like=None)
f, Welch = signal.welch(Txt, fs=1.0, window="hann", nperseg=None, noverlap=None, nfft=1024, detrend="constant", return_onesided=True, scaling="density", axis=-1, average="mean")
Colormap = plt.get_cmap("plasma")
Segment_Colormap = Colormap(np.linspace(0, 1, len(Subjects)))
plt.plot(f, Welch, c=Segment_Colormap[Subjects.index(Subject)], linewidth=2.5)
plt.show()
The following column is a shortened example of the input data for Subject1 (the .txt file input).
-0.315046 -0.0641514 -0.921697 -0.976516 -0.546322 -0.857073 -1.59524 -0.447158 0.559437 -0.0138681 0.308289 1.26514 1.4109 1.65135 0.70157 0.454801 1.20363 1.36459 0.26468 -0.487983 -0.181889 0.095632 -0.825288
