カラードノイズ生成ライブラリの解析

前提・実現したいこと

カラードノイズ生成プログラムの解析

発生している問題・エラーメッセージ

カラードノイズを使用したいと思っており、カラードノイズ生成ライブラリなるものを見つけたのですが、そのライブラリに載っていた例を実行した時のグラフの縦軸、横軸が何なのか(単位は何か)がわからず困っています。これがライブラリ載っていたサイトです。カラードノイズ生成ライブラリ
いかに示すのが実行時に出力されたグラフです。

該当のソースコード

これがまずカラードノイズ生成のライブラリです。

python
1"""Generate colored noise."""
2
3from numpy import sqrt, newaxis
4from numpy.fft import irfft, rfftfreq
5from numpy.random import normal
6from numpy import sum as npsum
7
8
9def powerlaw_psd_gaussian(exponent, size, fmin=0):
10    """Gaussian (1/f)**beta noise.
11
12    Based on the algorithm in:
13    Timmer, J. and Koenig, M.:
14    On generating power law noise.
15    Astron. Astrophys. 300, 707-710 (1995)
16
17    Normalised to unit variance
18
19    Parameters:
20    -----------
21
22    exponent : float
23        The power-spectrum of the generated noise is proportional to
24
25        S(f) = (1 / f)**beta
26        flicker / pink noise:   exponent beta = 1
27        brown noise:            exponent beta = 2
28
29        Furthermore, the autocorrelation decays proportional to lag**-gamma
30        with gamma = 1 - beta for 0 < beta < 1.
31        There may be finite-size issues for beta close to one.
32
33    shape : int or iterable
34        The output has the given shape, and the desired power spectrum in
35        the last coordinate. That is, the last dimension is taken as time,
36        and all other components are independent.
37
38    fmin : float, optional
39        Low-frequency cutoff.
40        Default: 0 corresponds to original paper. It is not actually
41        zero, but 1/samples.
42
43    Returns
44    -------
45    out : array
46        The samples.
47
48
49    Examples:
50    ---------
51
52    # generate 1/f noise == pink noise == flicker noise
53    >>> import colorednoise_output as cn
54    >>> y = cn.powerlaw_psd_gaussian(1, 5)
55    """
56    
57    # Make sure size is a list so we can iterate it and assign to it.
58    try:
59        size = list(size)
60    except TypeError:
61        size = [size]
62    
63    # The number of samples in each time series
64    samples = size[-1]
65    
66    # Calculate Frequencies (we asume a sample rate of one)
67    # Use fft functions for real output (-> hermitian spectrum)
68    f = rfftfreq(samples)
69    
70    # Build scaling factors for all frequencies
71    s_scale = f
72    fmin = max(fmin, 1./samples) # Low frequency cutoff
73    ix   = npsum(s_scale < fmin)   # Index of the cutoff
74    if ix and ix < len(s_scale):
75        s_scale[:ix] = s_scale[ix]
76    s_scale = s_scale**(-exponent/2.)
77    
78    # Calculate theoretical output standard deviation from scaling
79    w      = s_scale[1:].copy()
80    w[-1] *= (1 + (samples % 2)) / 2. # correct f = +-0.5
81    sigma = 2 * sqrt(npsum(w**2)) / samples
82    
83    # Adjust size to generate one Fourier component per frequency
84    size[-1] = len(f)
85
86    # Add empty dimension(s) to broadcast s_scale along last
87    # dimension of generated random power + phase (below)
88    dims_to_add = len(size) - 1
89    s_scale     = s_scale[(newaxis,) * dims_to_add + (Ellipsis,)]
90    
91    # Generate scaled random power + phase
92    sr = normal(scale=s_scale, size=size)
93    si = normal(scale=s_scale, size=size)
94    
95    # If the signal length is even, frequencies +/- 0.5 are equal
96    # so the coefficient must be real.
97    if not (samples % 2): si[...,-1] = 0
98    
99    # Regardless of signal length, the DC component must be real
100    si[...,0] = 0
101    
102    # Combine power + corrected phase to Fourier components
103    s  = sr + 1J * si
104    
105    # Transform to real time series & scale to unit variance
106    y = irfft(s, n=samples, axis=-1) / sigma
107    
108    return y
109

次にこれがサンプルとして乗っていた例です。
これを実行した際に出力されたグラフが上に示したものなのですが、グラフの横軸縦軸が、何を表すのかがわかっていない状態です。

python
1from matplotlib import mlab
2from matplotlib import pylab as plt
3import colorednoise as cn
4beta = 1 # the exponent
5samples = 2**18 # number of samples to generate
6y = cn.powerlaw_psd_gaussian(beta, samples)
7
8# optionally plot the Power Spectral Density with Matplotlib
9s, f = mlab.psd(y, NFFT=2**13)
10plt.loglog(f,s)
11plt.grid(True)
12plt.show()