我正在尝试使用 SciPy 的周期图和 AstroPy Lomb-Scargle 周期图来计算数据的周期图 - 周期图在任何地方都匹配,除了接近最小频率的频率,如我的图中所示。这些是数值模拟的结果。
根据观测数据,我预计在 0 附近有一个强信号。因此,SciPy 周期图结果看起来比 Lomb-Scargle 周期图在物理上更合理。
我还没弄清楚为什么以及如何使它们相似。任何见解都深表赞赏。
下面是重现我的绘图的代码。
来自标准 SciPy 周期图:
从 Lomb-Scargle 周期图:
from astropy.timeseries import LombScargle
import numpy as np
import pandas as pd
from scipy import signal
import requests
import matplotlib.pyplot as plt
def plot_periodogram(x,y,N_freq,min_freq,max_freq,height_threshold,periodogram_type):
fig, ax = plt.subplots(figsize=(12,8))
if periodogram_type == 'periodogram':
dx = np.mean(np.diff(x)) # Assume x is uniformly sampled
fs = 1 / dx
freq, power_periodogram = signal.periodogram(y,fs,scaling="spectrum",nfft=N_freq,
return_onesided=True,detrend='constant')
power_max = power_periodogram[~np.isnan(power_periodogram)].max()
plt.plot(freq, power_periodogram/power_max,linestyle="solid",color="black",linewidth=2)
filename = "PowerSpectrum"
else:
freq = np.linspace(min_freq,max_freq,N_freq)
ls= LombScargle(x, y,normalization='psd',nterms=1)
power_periodogram= ls.power(freq)
power_max = power_periodogram[~np.isnan(power_periodogram)].max()
false_alarm_probabilities = [0.01,0.05]
periodogram_peak_height= ls.false_alarm_level(false_alarm_probabilities,minimum_frequency=min_freq,
maximum_frequency=max_freq,method='bootstrap')
filename = "PowerSpectrum_LombScargle"
plt.plot(freq, power_periodogram/power_max,linestyle="solid",color="black",linewidth=2)
plt.axhline(y=periodogram_peak_height[0]/power_max, color='black', linestyle='--')
plt.axhline(y=periodogram_peak_height[1]/power_max, color='black', linestyle='-')
peaks_index, properties = signal.find_peaks(power_periodogram/power_max, height=height_threshold)
peak_values = properties['peak_heights']
peak_power_freq = freq[peaks_index]
for i in range(len(peak_power_freq)):
plt.axvline(x = peak_power_freq[i],color = 'red',linestyle='--')
ax.text(peak_power_freq[i]+0.05, 0.95, str(round(peak_power_freq[i],2)), color='red',ha='left', va='top', rotation=0,transform=ax.get_xaxis_transform())
fig.patch.set_alpha(1)
plt.ylabel('Spectral Power',fontsize=20)
plt.xlabel('Spatial Frequency', fontsize=20)
plt.grid(True)
plt.xlim(left=min_freq,right=max_freq)
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.savefig(filename,bbox_inches='tight')
plt.show()
# URL of the CSV file on Pastebin
url = 'https://pastebin.com/raw/uFi8WPvJ'
# Fetch the raw data from the URL
response = requests.get(url)
# Check if the request was successful
if response.status_code == 200:
# Decode the response content to text
data = response.text
# Save the data to a CSV file
with open('data.csv', 'w') as f:
f.write(data)
df =pd.read_csv('data.csv',sep=',',comment='%', names=['x', 'Bphi','r','theta'])
x = df['x'].values
y = df['Bphi'].values
# https://stackoverflow.com/questions/37540782/delete-nan-and-corresponding-elements-in-two-same-length-array
indices = np.logical_not(np.logical_or(np.isnan(x), np.isnan(y)))
x = x[indices]
y = y[indices]
y = y - np.mean(y)
N_freq = 10000
min_freq = 0.001;
max_freq = 4.0
height_threshold =0.7
plot_periodogram(x,y,N_freq,min_freq,max_freq,height_threshold,"periodogram")
plot_periodogram(x,y,N_freq,min_freq,max_freq,height_threshold,"ls")
我对代码进行了很多修改,并意识到我什至没有检查数据的样子:
数据不仅有偏移,还有趋势。在任何类型的频率变换之前都需要将其删除。
因此我使用了以下内容:
df =pd.read_csv('data.csv',sep=',',comment='%', names=['x', 'Bphi','r','theta'])
df.dropna(inplace = True)
x = df['x'].values
y = signal.detrend(df['Bphi'].values)
结果并不相同,但非常相似。
Scipy 方法:
天向性方法:
我建议深入研究这两个函数的文档。完善此方法结束后,您可以使用
np.allclose()
检查结果是否可接受。