<p>This is a partial answer based on the discussion in comments below the question.</p>
<p>If one already has an &quot;amplitude spectrum&quot; and one wants to convert to a power spectrum, all you have to do is take the absolute value of the amplitude and square it. In Python (introduced in the question) that's just <code>np.abs(amplitude)**2</code>.</p>
<p>When you take the Discrete Fourier transform of a time series it is returned as a complex array, you need both real and imaginary components to retrieve both the absolute value of the amplitude spectrum and its phase spectrum.</p>
<p>Sometimes you can see &quot;amplitude&quot; used for the absolute value of amplitude, but it doesn't usually hurt if you take the absolute value again to be sure.</p>
<p>I've combined the signals of random noise with <span class="math-container">$\sigma$</span> = 1, a sinusoidal 50 Hz electronic interference noise signal, and a slightly phase-jittery square wave with an amplitude of 0.5 and a period of 12 milliseconds or a frequency of about 83 Hz.</p>
<p>I've kept the average near zero so we don't get distracted by a giant peak at 0 Hz, and I've shifted zero frequency to the center of the spectrum using <code>np.fft.fftshift</code> so we don't forget that negative frequencies exist in Fourier space.</p>
<p>We can see that the FT has both real and imaginary parts and squaring the absolute value gives us power.</p>
<p>For the time series plots I only show the first 200 milliseconds for clarity, the series has 1,000 time intervals of 1 millisecond each.</p>
<p>I've also used <code>np.fft.fftfreq</code> to automatically generate the frequency assigned to each bin in the FT.</p>
<hr />
<p><a href="https://i.sstatic.net/5flYE.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/5flYE.png" alt="enter image description here" /></a></p>
<pre><code>import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42) # this way you can compare your output to mine directly

N = 2000
nt = np.arange(N)
dt = 0.001 # one millisecond per bin
time = nt * dt

twopi = 2 * np.pi

noise = np.random.randn(N) # normal distribution with sigma = 1
frequency_1 = 50
signal_1 = 0.2 * np.cos(twopi * time * frequency_1)

# jitter = (7 * np.random.random(N)).astype(int) - 3
interval_2 = 12
q_int = int(interval_2 / 4)
ints = np.arange(0, N, interval_2)
jittered = np.random.randint(-2, 2, len(ints)) + ints
signal_2 = np.zeros_like(nt, dtype='float')
for j in jittered:
    signal_2[j-q_int:j+q_int+1] += 1 # I suppose one could convolute or sum over np.shift
signal_2 -= signal_2.mean()

combined = noise + signal_1 + signal_2

ft = np.fft.fft(combined)
freqs = np.fft.fftfreq(N, dt)

ft_shifted = np.fft.fftshift(ft)
power_shifted = np.abs(ft_shifted)**2
freq_shifted = np.fft.fftshift(freqs)

things, names = ((noise, signal_1, signal_2, combined),
                 ('noise', 'signal_1', 'signal_2', 'combined'))
m = len(things)
fig, axes = plt.subplots(m+2)
for i, (thing, name, ax) in enumerate(zip(things, names, axes[:-2])):
    ax.plot(time, thing)
    ax.plot(time, np.zeros_like(time), '-k', linewidth=0.5)
    ax.set_xlim(0, 0.2)
    ax.set_ylabel(name)
axes[-2].plot(freq_shifted, ft_shifted.real)
axes[-2].plot(freq_shifted, ft_shifted.imag)
axes[-1].plot(freq_shifted, power_shifted)
axes[-1].set_xlabel('frequency (Hz)')
for ax in axes[-2:]:
    ax.set_xlim(-150, 150)
plt.show()
</code></pre>

CC BY-SA 4.0
2021-03-27T10:30:17.267Z
I did some digging with your code and manage to get the answer altho. the 'tick markers' are in the negative side (I'll try to fix it). Binned instead to dt = 0.1. Thanks by a lot.
CGHA
CC BY-SA 4.0
2021-03-27T11:42:27.927Z
I can see your deleted answer (nothing is really deleted in Stack Exchange) and I don't know why somebody down-voted it. I think we have a user who is relatively new to Astronomy SE who is down and close-voting more actively than most folks. *It is always okay* to answer your own question, and my answer self-identifies as only a "partial answer".  Why not undelete your answer, I'll give it an up vote and you can just edit the first sentence to sound more like you are posting an answer for others to read instead of a reply to me? Just explain a bit why it addresses the question.
uhoh
CC BY-SA 4.0
2021-03-27T04:30:24.200Z
2021-03-27T04:30:24.200Z
uhoh
<p>So after further digging and following the partial answer above, I came up with this code/answer which should suffice if you got the 'amplitude spectrum' right away and works.</p>
<pre><code>import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

df = pd.read_csv('FLSP.csv', index_col = False) 
with pd.option_context('display.float_format', '{:0.20f}'.format): #optional
    print(df) #Data Frame

N = df.shape[0]
dt = 0.1 #binning at one second

ft = np.fft.fft(data)
freqs = np.fft.fftfreq(N, dt)

ft_shifted = np.fft.fftshift(ft)
power_shifted = np.abs(data['Amplitude'])**2
freq_shifted = np.fft.fftshift(freqs)

plt.plot(freq_shifted, power_shifted, label=&quot;Power Spectrum&quot;) 
plt.xlabel('Frequency [Hz]')
plt.ylabel('Power')
</code></pre>
<p>Which should show something like this;</p>
<p><a href="https://i.sstatic.net/GuHBB.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/GuHBB.png" alt="enter image description here" /></a></p>
<p>It's close yet I know there's still a lot of work to do. Unfortunately, the data is in the &quot;negative side&quot; for some reason (hence with the presented tick markers) which I'll address later on.</p>
<p>Definitely, thanks by a lot c:</p>

CC BY-SA 4.0
2021-03-27T09:01:25.600Z
2021-03-27T12:36:04.880Z
2021-03-27T12:36:04.880Z
CGHA
<p>So I already did a frequency analysis of the SRv-type's light curve  using <a href="http://www.period04.net/" rel="nofollow noreferrer">Period04</a>. I already have the 'periodicities' of the object/star of choice which is quite hard since the latter/other periods are tucked away by yearly aliases <em>(thank goodness that's a rough week)</em>.</p>
<p>The problem, however, is that the mentioned-software doesn't allow the user to acquire the Power Spectrum but only the Fourier data (as of the moment) which is important for visualizing the non-coherent signals. May I kindly ask if it's possible to deduce/extract the power spectrum from the fourier data using Lomb-Scargle periodogram? And if there's an existing python tool/code to it?</p>
<p>The current code, which is below bog standard I'm afraid, just visualizes the Fourier spectrum;</p>
<pre><code>import pandas as pd
import matplotlib.pyplot as plt
with pd.option_context('display.float_format', '{:0.20f}'.format):
    data = pd.read_csv('FLSP.csv')
    print(data) #shows Frequencies/Amplitudes

x = data['Frequency']
y = data['Amplitude']

plt.xlabel('Frequency')
plt.ylabel('Amplitude')

plt.plot(x, y)
</code></pre>
<p>Which should give you something (as useless) like this;</p>
<p><a href="https://i.sstatic.net/Sf0d8.jpg" rel="nofollow noreferrer"><img src="https://i.sstatic.net/Sf0d8.jpg" alt="enter image description here" /></a></p>
<p>It's the Long Secondary Period (LSP) of the variable star at ~2350 days. Attach here is the <a href="https://drive.google.com/file/d/1kIEDG518hW8thIXCtLkr-UJMJOT-jTNm/view" rel="nofollow noreferrer">fourier data</a> if you want to check out.</p>
<p>Help at this point is sincerely appreciated. Thnx a lot and clear skies.</p>
<p>Note: I tried using Astropy's and Scipy's Lomb-Scargle, and their implementations require the observations; light curve, so I'm stuck on 'looking' at the fourier spectrum...suggesting to me that's its back to square one (hoping not lol).</p>

CC BY-SA 4.0
2021-03-27T01:42:30.433Z
I understand wanting to post the same question on multiple sites, we all have that experience. But cross-posting is strongly discouraged for reasons including *answer fragmentation*. Future readers may come across one set of answers and not always find the other set. So we really should choose the best one site to post, and delete other copies. That said, I wonder if the best one site might be [Signal Processing SE](https://dsp.stackexchange.com/) rather than Astronomy SE or Stack Oveflow?
uhoh
CC BY-SA 4.0
2021-03-27T02:23:15.433Z
Thanks for the response and my bad. I just thought it was astronomy related since we're talking about light curve analysis of stars not just ordinary signals, primarily about period/frequency analysis, which is done by multiple astro peeps and hoping a few response from them...if they're ever here at lol.
CGHA
CC BY-SA 4.0
2021-03-27T02:46:31.513Z
Ya I agree, I think that this question is just fine here, on-topic and may well receive a a good answer soon. If it were me, I'd delete the Stack Overflow copy at some point, at the latest when an answer is posted here.
uhoh
CC BY-SA 4.0
2021-03-27T02:49:29.783Z
Are you just asking how to change an amplitude spectrum ("Fourier data") into a power spectrum? Isn't that what `np.abs(amplitude)**2' does? If that's what it is that you want to do, then I'm not sure what Lomb-Scargle has to do with that. Is it possible that you are making the problem more complicated than it needs to be?
uhoh
CC BY-SA 4.0
2021-03-27T02:55:43.767Z
Sincere there's multiple questions, I'll break it down into 3 answers: 
1. Yes that's pretty much what I wanted to be -> change/convert the amplitude (fourier) spectrum to power spectrum.
2. I would check that np.abs(amplitude)**2 thing but isn't that a part of DFT/FFT instead of L-S periodogram? 
3. Yes it's possible that I'm making problem more worse than it suppose to be haha. Unfortunately, being not well-verse in Python like another perpetrator to that.
CGHA
CC BY-SA 4.0
2021-03-27T03:32:43.313Z
Don't get me wrong (even tho. I'm not skilled at python as anyone), I'm also considering FT functions, primarily the NFFT (or the non-equal 'version' of FFT). But, the reason why I'm trying to incorporate L-S periodogram for this one, aside from being the most utilized statistical method in light curve data, is the capacity to maximize the likelihood of peak detection and residual detection (since the object 'pulses' for 3 times at different scales; 2 of the detection are tucked away by the aliases) and noise is distributed relatively similar to a chi-squared manner.
CGHA
CC BY-SA 4.0
2021-03-27T04:34:26.623Z
I've posted a partial answer to get things started. Use `import numpy as np` Numpy or Numerical Python is an incredibly useful thing to know, its used everywhere and inside most Python packages that you use. If you are going to use math, then get to know Numpy, even if you stick to standard packages for serious work (it's so easy to make mistakes ourselves) it's a good thing to know.  *All that being said* you may still need an another answer if you want to convert amplitude to power. I think you just square the absolute value and square the units, but let's see what others think.
uhoh
<data-analysis><variable-star>
Frequency analysis on Semiperiodic object(s) using Lomb-Scargle
CGHA
# Frequency analysis on Semiperiodic object(s) using Lomb-Scargle

(Asked by: username_0 on 2021-03-27)

So I already did a frequency analysis of the SRv-type's light curve using [Period04](http://www.period04.net/). I already have the 'periodicities' of the object/star of choice which is quite hard since the latter/other periods are tucked away by yearly aliases *(thank goodness that's a rough week)*.

The problem, however, is that the mentioned-software doesn't allow the user to acquire the Power Spectrum but only the Fourier data (as of the moment) which is important for visualizing the non-coherent signals. May I kindly ask if it's possible to deduce/extract the power spectrum from the fourier data using Lomb-Scargle periodogram? And if there's an existing python tool/code to it?

The current code, which is below bog standard I'm afraid, just visualizes the Fourier spectrum;

```
import pandas as pd
import matplotlib.pyplot as plt
with pd.option_context('display.float_format', '{:0.20f}'.format):
    data = pd.read_csv('FLSP.csv')
    print(data) #shows Frequencies/Amplitudes

x = data['Frequency']
y = data['Amplitude']

plt.xlabel('Frequency')
plt.ylabel('Amplitude')

plt.plot(x, y)

```

Which should give you something (as useless) like this;

[![enter image description here](https://i.sstatic.net/Sf0d8.jpg)](https://i.sstatic.net/Sf0d8.jpg)

It's the Long Secondary Period (LSP) of the variable star at ~2350 days. Attach here is the [fourier data](https://drive.google.com/file/d/1kIEDG518hW8thIXCtLkr-UJMJOT-jTNm/view) if you want to check out.

Help at this point is sincerely appreciated. Thnx a lot and clear skies.

Note: I tried using Astropy's and Scipy's Lomb-Scargle, and their implementations require the observations; light curve, so I'm stuck on 'looking' at the fourier spectrum...suggesting to me that's its back to square one (hoping not lol).

**Tags:** data-analysis, variable-star

> **Comments:**
> **username_1**: I understand wanting to post the same question on multiple sites, we all have that experience. But cross-posting is strongly discouraged for reasons including *answer fragmentation*. Future readers may come across one set of answers and not always find the other set. So we really should choose the best one site to post, and delete other copies. That said, I wonder if the best one site might be [Signal Processing SE](https://dsp.stackexchange.com/) rather than Astronomy SE or Stack Oveflow?  
> **username_0**: Thanks for the response and my bad. I just thought it was astronomy related since we're talking about light curve analysis of stars not just ordinary signals, primarily about period/frequency analysis, which is done by multiple astro peeps and hoping a few response from them...if they're ever here at lol.  
> **username_1**: Ya I agree, I think that this question is just fine here, on-topic and may well receive a a good answer soon. If it were me, I'd delete the Stack Overflow copy at some point, at the latest when an answer is posted here.  
> **username_1**: Are you just asking how to change an amplitude spectrum ("Fourier data") into a power spectrum? Isn't that what `np.abs(amplitude)**2' does? If that's what it is that you want to do, then I'm not sure what Lomb-Scargle has to do with that. Is it possible that you are making the problem more complicated than it needs to be?  
> **username_0**: Sincere there's multiple questions, I'll break it down into 3 answers:  1. Yes that's pretty much what I wanted to be -> change/convert the amplitude (fourier) spectrum to power spectrum. 2. I would check that np.abs(amplitude)**2 thing but isn't that a part of DFT/FFT instead of L-S periodogram?  3. Yes it's possible that I'm making problem more worse than it suppose to be haha. Unfortunately, being not well-verse in Python like another perpetrator to that.  
> **username_0**: Don't get me wrong (even tho. I'm not skilled at python as anyone), I'm also considering FT functions, primarily the NFFT (or the non-equal 'version' of FFT). But, the reason why I'm trying to incorporate L-S periodogram for this one, aside from being the most utilized statistical method in light curve data, is the capacity to maximize the likelihood of peak detection and residual detection (since the object 'pulses' for 3 times at different scales; 2 of the detection are tucked away by the aliases) and noise is distributed relatively similar to a chi-squared manner.  
> **username_1**: I've posted a partial answer to get things started. Use `import numpy as np` Numpy or Numerical Python is an incredibly useful thing to know, its used everywhere and inside most Python packages that you use. If you are going to use math, then get to know Numpy, even if you stick to standard packages for serious work (it's so easy to make mistakes ourselves) it's a good thing to know.  *All that being said* you may still need an another answer if you want to convert amplitude to power. I think you just square the absolute value and square the units, but let's see what others think.  

---

## 2 answers

### Answer by username_1 on 2021-03-27. Score: 1 | Accepted answer

This is a partial answer based on the discussion in comments below the question.

If one already has an "amplitude spectrum" and one wants to convert to a power spectrum, all you have to do is take the absolute value of the amplitude and square it. In Python (introduced in the question) that's just `np.abs(amplitude)**2`.

When you take the Discrete Fourier transform of a time series it is returned as a complex array, you need both real and imaginary components to retrieve both the absolute value of the amplitude spectrum and its phase spectrum.

Sometimes you can see "amplitude" used for the absolute value of amplitude, but it doesn't usually hurt if you take the absolute value again to be sure.

I've combined the signals of random noise with $\sigma$ = 1, a sinusoidal 50 Hz electronic interference noise signal, and a slightly phase-jittery square wave with an amplitude of 0.5 and a period of 12 milliseconds or a frequency of about 83 Hz.

I've kept the average near zero so we don't get distracted by a giant peak at 0 Hz, and I've shifted zero frequency to the center of the spectrum using `np.fft.fftshift` so we don't forget that negative frequencies exist in Fourier space.

We can see that the FT has both real and imaginary parts and squaring the absolute value gives us power.

For the time series plots I only show the first 200 milliseconds for clarity, the series has 1,000 time intervals of 1 millisecond each.

I've also used `np.fft.fftfreq` to automatically generate the frequency assigned to each bin in the FT.

---

[![enter image description here](https://i.sstatic.net/5flYE.png)](https://i.sstatic.net/5flYE.png)

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

np.random.seed(42) # this way you can compare your output to mine directly

N = 2000
nt = np.arange(N)
dt = 0.001 # one millisecond per bin
time = nt * dt

twopi = 2 * np.pi

noise = np.random.randn(N) # normal distribution with sigma = 1
frequency_1 = 50
signal_1 = 0.2 * np.cos(twopi * time * frequency_1)

# jitter = (7 * np.random.random(N)).astype(int) - 3
interval_2 = 12
q_int = int(interval_2 / 4)
ints = np.arange(0, N, interval_2)
jittered = np.random.randint(-2, 2, len(ints)) + ints
signal_2 = np.zeros_like(nt, dtype='float')
for j in jittered:
    signal_2[j-q_int:j+q_int+1] += 1 # I suppose one could convolute or sum over np.shift
signal_2 -= signal_2.mean()

combined = noise + signal_1 + signal_2

ft = np.fft.fft(combined)
freqs = np.fft.fftfreq(N, dt)

ft_shifted = np.fft.fftshift(ft)
power_shifted = np.abs(ft_shifted)**2
freq_shifted = np.fft.fftshift(freqs)

things, names = ((noise, signal_1, signal_2, combined),
                 ('noise', 'signal_1', 'signal_2', 'combined'))
m = len(things)
fig, axes = plt.subplots(m+2)
for i, (thing, name, ax) in enumerate(zip(things, names, axes[:-2])):
    ax.plot(time, thing)
    ax.plot(time, np.zeros_like(time), '-k', linewidth=0.5)
    ax.set_xlim(0, 0.2)
    ax.set_ylabel(name)
axes[-2].plot(freq_shifted, ft_shifted.real)
axes[-2].plot(freq_shifted, ft_shifted.imag)
axes[-1].plot(freq_shifted, power_shifted)
axes[-1].set_xlabel('frequency (Hz)')
for ax in axes[-2:]:
    ax.set_xlim(-150, 150)
plt.show()

```

> **Comments:**
> **username_0**: I did some digging with your code and manage to get the answer altho. the 'tick markers' are in the negative side (I'll try to fix it). Binned instead to dt = 0.1. Thanks by a lot.  
> **username_1**: I can see your deleted answer (nothing is really deleted in Stack Exchange) and I don't know why somebody down-voted it. I think we have a user who is relatively new to Astronomy SE who is down and close-voting more actively than most folks. *It is always okay* to answer your own question, and my answer self-identifies as only a "partial answer".  Why not undelete your answer, I'll give it an up vote and you can just edit the first sentence to sound more like you are posting an answer for others to read instead of a reply to me? Just explain a bit why it addresses the question.  

---

### Answer by username_0 on 2021-03-27. Score: 0

So after further digging and following the partial answer above, I came up with this code/answer which should suffice if you got the 'amplitude spectrum' right away and works.

```
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

df = pd.read_csv('FLSP.csv', index_col = False) 
with pd.option_context('display.float_format', '{:0.20f}'.format): #optional
    print(df) #Data Frame

N = df.shape[0]
dt = 0.1 #binning at one second

ft = np.fft.fft(data)
freqs = np.fft.fftfreq(N, dt)

ft_shifted = np.fft.fftshift(ft)
power_shifted = np.abs(data['Amplitude'])**2
freq_shifted = np.fft.fftshift(freqs)

plt.plot(freq_shifted, power_shifted, label="Power Spectrum") 
plt.xlabel('Frequency [Hz]')
plt.ylabel('Power')

```

Which should show something like this;

[![enter image description here](https://i.sstatic.net/GuHBB.png)](https://i.sstatic.net/GuHBB.png)

It's close yet I know there's still a lot of work to do. Unfortunately, the data is in the "negative side" for some reason (hence with the presented tick markers) which I'll address later on.

Definitely, thanks by a lot c:


---

