0
$\begingroup$

So I already did a frequency analysis of the SRv-type's light curve using Period04. 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

It's the Long Secondary Period (LSP) of the variable star at ~2350 days. Attach here is the fourier data 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).

$\endgroup$
7
  • 1
    $\begingroup$ 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 rather than Astronomy SE or Stack Oveflow? $\endgroup$ – uhoh Mar 27 at 1:42
  • 1
    $\begingroup$ 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. $\endgroup$ – CGHA Mar 27 at 2:23
  • $\begingroup$ 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. $\endgroup$ – uhoh Mar 27 at 2:46
  • $\begingroup$ 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? $\endgroup$ – uhoh Mar 27 at 2:49
  • 1
    $\begingroup$ 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. $\endgroup$ – CGHA Mar 27 at 2:55
1
$\begingroup$

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

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()
$\endgroup$
2
  • 1
    $\begingroup$ 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. $\endgroup$ – CGHA Mar 27 at 10:30
  • 1
    $\begingroup$ 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. $\endgroup$ – uhoh Mar 27 at 11:42
0
$\begingroup$

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

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:

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.