Несоответствие между периодограммой, рассчитанной с помощью периодограммы SciPy, и периодограммы AstroPy Lomb Scargle на низких частотах

Я пытаюсь вычислить периодограмму своих данных, используя как периодограмму SciPy, так и периодограмму AstroPy Периодограмму Ломба-Скаргла — периодограмма совпадает везде, кроме частот, близких к минимальной частоте, как показано на моих графиках. Это результаты численного моделирования.

Основываясь на данных наблюдений, я ожидаю сильный сигнал около 0. Следовательно, результаты периодограммы SciPy выглядят более физически правдоподобными, чем периодограмма Ломба-Скаргла.

Я так и не понял, почему и как сделать их похожими. Любое понимание глубоко ценится.

Ниже приведен код для воспроизведения моих графиков.

Из стандартной периодограммы SciPy: Из периодограммы Ломб-Скаргла:

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")

   

это значительное смещение постоянного тока. Вы проверяли, удалил ли метод Ломб-Скаргла среднее значение из данных?

Tino D 02.06.2024 09:55

из scipy документации. По умолчанию для удаления тренда установлено значение «константа», которое удаляет среднее значение из данных. либо укажите его как ложное, либо удалите среднее значение самостоятельно перед Ломбом-Скарглом...

Tino D 02.06.2024 09:59

@TinoD Я вручную удаляю среднее значение из данных перед вызовом любой функции. Вы найдете его в конце моего фрагмента кода y = y - np.mean(y)

Prav001 02.06.2024 10:00
Почему в Python есть оператор "pass"?
Почему в Python есть оператор "pass"?
Оператор pass в Python - это простая концепция, которую могут быстро освоить даже новички без опыта программирования.
Некоторые методы, о которых вы не знали, что они существуют в Python
Некоторые методы, о которых вы не знали, что они существуют в Python
Python - самый известный и самый простой в изучении язык в наши дни. Имея широкий спектр применения в области машинного обучения, Data Science,...
Основы Python Часть I
Основы Python Часть I
Вы когда-нибудь задумывались, почему в программах на Python вы видите приведенный ниже код?
LeetCode - 1579. Удаление максимального числа ребер для сохранения полной проходимости графа
LeetCode - 1579. Удаление максимального числа ребер для сохранения полной проходимости графа
Алиса и Боб имеют неориентированный граф из n узлов и трех типов ребер:
Оптимизация кода с помощью тернарного оператора Python
Оптимизация кода с помощью тернарного оператора Python
И последнее, что мы хотели бы показать вам, прежде чем двигаться дальше, это
Советы по эффективной веб-разработке с помощью Python
Советы по эффективной веб-разработке с помощью Python
Как веб-разработчик, Python может стать мощным инструментом для создания эффективных и масштабируемых веб-приложений.
1
3
87
1
Перейти к ответу Данный вопрос помечен как решенный

Ответы 1

Ответ принят как подходящий

Я много возился с кодом и понял, что даже не проверял, как выглядят данные:

Данные имеют не только смещение, но и тенденцию. Это необходимо удалить перед любым преобразованием частоты.

Поэтому я использовал следующее:

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)

Результаты не идентичны, но очень похожи.

Сциповый подход:

Астропический подход:

Я бы порекомендовал очень глубоко изучить документацию обеих функций. В конце доработки этого подхода вы можете использовать np.allclose(), чтобы проверить, приемлемы ли результаты.

Снятие тренда в scipy по умолчанию равно linear, а в периодограмме в качестве удаления тренда по умолчанию используется constant. Вероятно, из-за этого небольшие различия в сюжете.

dankal444 02.06.2024 14:25

@dankal444 может быть, я думаю, что очень глубокое чтение документации может помочь в этой ситуации.

Tino D 02.06.2024 14:45

@TinoD Спасибо! Я пытался исключить тренд из данных, но использовал type = 'constant' вместо значения по умолчанию type='linear'. Когда я использую первый вариант, я получаю те же результаты, что и в исходном сюжете. Не знаю, почему это происходит, но думаю, мне нужно углубиться в документацию.

Prav001 02.06.2024 22:04

Другие вопросы по теме