數(shù)據(jù)分析：時間序列分析：時間序列的譜分析

數(shù)據(jù)分析：時間序列分析：時間序列的譜分析1數(shù)據(jù)分析：時間序列分析：時間序列的譜分析1.1時間序列分析簡介1.1.1時間序列數(shù)據(jù)的特點(diǎn)時間序列數(shù)據(jù)是指在一系列時間點(diǎn)上收集的數(shù)據(jù)點(diǎn)，這些數(shù)據(jù)點(diǎn)通常按照時間順序排列，反映了某個現(xiàn)象或過程隨時間變化的趨勢。時間序列數(shù)據(jù)具有以下特點(diǎn)：順序性：數(shù)據(jù)點(diǎn)的順序至關(guān)重要，因?yàn)樗鼈兎从沉藭r間的流逝。周期性：許多時間序列數(shù)據(jù)展現(xiàn)出周期性的模式，如季節(jié)性波動。趨勢性：數(shù)據(jù)可能隨時間呈現(xiàn)上升、下降或平穩(wěn)的趨勢。隨機(jī)性：除了明顯的趨勢和周期，時間序列數(shù)據(jù)還可能包含隨機(jī)波動。自相關(guān)性：當(dāng)前時間點(diǎn)的數(shù)據(jù)值與過去時間點(diǎn)的數(shù)據(jù)值之間可能存在相關(guān)性。1.1.2時間序列分析的基本步驟時間序列分析的基本步驟包括：數(shù)據(jù)預(yù)處理：清洗數(shù)據(jù)，處理缺失值，平滑數(shù)據(jù)等。探索性分析：可視化數(shù)據(jù)，檢查趨勢、季節(jié)性和周期性。模型選擇：根據(jù)數(shù)據(jù)特性選擇合適的模型，如ARIMA、狀態(tài)空間模型等。參數(shù)估計(jì)：使用數(shù)據(jù)估計(jì)模型參數(shù)。模型診斷：檢查模型的殘差是否滿足模型假設(shè)。預(yù)測：使用模型對未來數(shù)據(jù)進(jìn)行預(yù)測。模型驗(yàn)證：通過回測或使用測試集數(shù)據(jù)驗(yàn)證模型的預(yù)測能力。示例：使用Python進(jìn)行時間序列分析#導(dǎo)入必要的庫

importpandasaspd

importnumpyasnp

importmatplotlib.pyplotasplt

fromstatsmodels.tsa.stattoolsimportadfuller

fromstatsmodels.tsa.seasonalimportseasonal_decompose

fromstatsmodels.tsa.arima.modelimportARIMA

#加載數(shù)據(jù)

data=pd.read_csv('sales_data.csv',index_col='Date',parse_dates=True)

#數(shù)據(jù)預(yù)處理

data=data.fillna(method='ffill')#使用前向填充處理缺失值

#探索性分析

plt.figure(figsize=(12,6))

plt.plot(data['Sales'])

plt.title('時間序列數(shù)據(jù)')

plt.xlabel('時間')

plt.ylabel('銷售額')

plt.show()

#檢查平穩(wěn)性

result=adfuller(data['Sales'])

print(f'ADFStatistic:{result[0]}')

print(f'p-value:{result[1]}')

#分解時間序列

decomposition=seasonal_decompose(data['Sales'],model='additive',period=12)

decomposition.plot()

plt.show()

#模型選擇與參數(shù)估計(jì)

model=ARIMA(data['Sales'],order=(1,1,1))

model_fit=model.fit()

#模型診斷

residuals=pd.DataFrame(model_fit.resid)

residuals.plot()

plt.show()

#預(yù)測

forecast=model_fit.forecast(steps=10)

print(forecast)

#模型驗(yàn)證

#使用回測數(shù)據(jù)驗(yàn)證模型預(yù)測能力在這個例子中，我們使用了Python的pandas庫來加載和預(yù)處理數(shù)據(jù)，matplotlib庫進(jìn)行數(shù)據(jù)可視化，statsmodels庫進(jìn)行時間序列分析，包括平穩(wěn)性檢驗(yàn)、季節(jié)性分解和ARIMA模型的構(gòu)建與預(yù)測。數(shù)據(jù)樣例假設(shè)sales_data.csv文件包含以下數(shù)據(jù)：Date,Sales

2020-01-01,100

2020-01-02,105

2020-01-03,110

...

2021-12-31,150這個數(shù)據(jù)集記錄了每天的銷售額，我們可以使用上述代碼進(jìn)行時間序列分析，包括數(shù)據(jù)預(yù)處理、探索性分析、模型選擇、參數(shù)估計(jì)、模型診斷、預(yù)測和模型驗(yàn)證。2數(shù)據(jù)分析：時間序列分析：譜分析基礎(chǔ)2.1傅里葉變換原理傅里葉變換是一種將時間序列信號轉(zhuǎn)換為頻率域表示的數(shù)學(xué)工具。它基于一個基本假設(shè)：任何復(fù)雜的信號都可以表示為不同頻率的正弦波和余弦波的線性組合。在時間序列分析中，傅里葉變換幫助我們識別信號中的周期性模式和頻率成分。2.1.1離散傅里葉變換(DFT)離散傅里葉變換是傅里葉變換在離散時間信號上的應(yīng)用。對于一個長度為N的離散時間序列xnX其中Xk是信號在頻率k上的頻譜分量，e2.1.2快速傅里葉變換(FFT)快速傅里葉變換是一種高效的算法，用于計(jì)算DFT。它通過將DFT分解為更小的DFT計(jì)算來減少計(jì)算量，從而大大提高了計(jì)算速度。在實(shí)際應(yīng)用中，F(xiàn)FT是進(jìn)行頻譜分析的首選方法。示例代碼importnumpyasnp

importmatplotlib.pyplotasplt

#創(chuàng)建一個時間序列數(shù)據(jù)

N=600

T=1.0/800.0

x=np.linspace(0.0,N*T,N)

y=np.sin(50.0*2.0*np.pi*x)+0.5*np.sin(80.0*2.0*np.pi*x)

#使用FFT進(jìn)行變換

yf=np.fft.fft(y)

xf=np.fft.fftfreq(N,T)

#繪制頻譜圖

plt.plot(xf,np.abs(yf))

plt.grid()

plt.xlim(0,100)

plt.xlabel('頻率')

plt.ylabel('幅度')

plt.show()這段代碼首先生成了一個包含兩個不同頻率的正弦波的時間序列數(shù)據(jù)。然后，使用numpy庫中的fft函數(shù)對數(shù)據(jù)進(jìn)行快速傅里葉變換，得到頻譜分量。最后，使用matplotlib庫繪制頻譜圖，顯示了信號中不同頻率的幅度。2.2功率譜密度概念功率譜密度（PowerSpectralDensity，PSD）是描述信號在不同頻率上的功率分布的統(tǒng)計(jì)量。它提供了信號能量在頻率域上的分布情況，對于分析信號的頻率特性非常有用。2.2.1計(jì)算PSDPSD可以通過計(jì)算信號的自相關(guān)函數(shù)的傅里葉變換得到，也可以直接通過對信號的傅里葉變換結(jié)果進(jìn)行平方并歸一化得到。在實(shí)際應(yīng)用中，通常使用后一種方法，因?yàn)樗苯忧矣?jì)算效率更高。示例代碼importnumpyasnp

importmatplotlib.pyplotasplt

fromscipy.signalimportwelch

#創(chuàng)建一個時間序列數(shù)據(jù)

N=1000

T=1.0/100.0

x=np.linspace(0.0,N*T,N)

y=np.sin(2.0*np.pi*10.0*x)+0.5*np.sin(2.0*np.pi*20.0*x)

#使用Welch方法計(jì)算PSD

frequencies,psd=welch(y,fs=100)

#繪制PSD圖

plt.semilogy(frequencies,psd)

plt.xlabel('頻率')

plt.ylabel('PSD')

plt.grid()

plt.show()在這個例子中，我們使用了scipy庫中的welch函數(shù)來計(jì)算PSD。welch函數(shù)使用了Welch方法，這是一種通過將信號分割成多個重疊的段，然后對每個段進(jìn)行FFT并平均結(jié)果來估計(jì)PSD的方法。這種方法可以減少估計(jì)的方差，提高PSD估計(jì)的準(zhǔn)確性。2.2.2解釋PSD圖PSD圖通常以對數(shù)刻度繪制，因?yàn)樾盘柕墓β史植纪缭蕉鄠€數(shù)量級。在上面的示例中，我們可以看到兩個峰值，分別對應(yīng)于信號中的10Hz和20Hz頻率成分。PSD圖的形狀和峰值位置可以幫助我們識別信號中的主要頻率成分，這對于診斷機(jī)械系統(tǒng)的振動、分析心電圖信號等應(yīng)用非常有用。通過以上內(nèi)容，我們了解了傅里葉變換原理和功率譜密度的概念，以及如何使用Python中的numpy和scipy庫來計(jì)算和可視化這些頻譜特性。這些工具和方法是進(jìn)行時間序列譜分析的基礎(chǔ)，可以幫助我們深入理解信號的頻率特性。3時間序列的譜分析方法3.1周期圖法周期圖法是時間序列譜分析中最直觀的方法之一，它通過計(jì)算時間序列的傅里葉變換來揭示序列中可能存在的周期性成分。周期圖（SpectralDensityPlot）顯示了不同頻率下的功率分布，幫助我們識別出序列中的主要頻率。3.1.1原理時間序列的周期圖是其自相關(guān)函數(shù)的傅里葉變換。對于一個時間序列xt，其周期圖II其中，γt3.1.2示例代碼假設(shè)我們有一個時間序列數(shù)據(jù)，我們使用Python的matplotlib和scipy庫來繪制其周期圖。importnumpyasnp

importmatplotlib.pyplotasplt

fromscipy.signalimportperiodogram

#生成一個具有周期性的模擬時間序列

np.random.seed(42)

t=np.linspace(0,10,1000,endpoint=False)

x=np.sin(2*np.pi*5*t)+np.random.normal(0,0.5,t.shape)

#計(jì)算周期圖

frequencies,spectrum=periodogram(x,fs=100)

#繪制周期圖

plt.figure(figsize=(10,5))

plt.semilogy(frequencies,spectrum)

plt.title('周期圖')

plt.xlabel('頻率')

plt.ylabel('功率譜密度')

plt.grid(True)

plt.show()3.1.3解釋在上述代碼中，我們首先生成了一個具有5Hz頻率的正弦波，并添加了一些隨機(jī)噪聲。然后，我們使用scipy.signal.periodogram函數(shù)計(jì)算了序列的周期圖。最后，我們使用matplotlib庫繪制了周期圖，其中頻率以Hz為單位，功率譜密度以對數(shù)刻度顯示。3.2自相關(guān)函數(shù)譜分析自相關(guān)函數(shù)譜分析是另一種常用的時間序列譜分析方法，它基于時間序列的自相關(guān)函數(shù)（AutocorrelationFunction,ACF）。自相關(guān)函數(shù)描述了序列中不同時間點(diǎn)上的值之間的相關(guān)性，而自相關(guān)函數(shù)的譜分析則揭示了這些相關(guān)性在頻率域中的分布。3.2.1原理自相關(guān)函數(shù)譜分析通過將自相關(guān)函數(shù)ρk進(jìn)行傅里葉變換來得到譜密度SS其中，ρk3.2.2示例代碼我們繼續(xù)使用Python的matplotlib和statsmodels庫來展示如何進(jìn)行自相關(guān)函數(shù)譜分析。importnumpyasnp

importmatplotlib.pyplotasplt

importstatsmodels.apiassm

#使用相同的模擬時間序列

t=np.linspace(0,10,1000,endpoint=False)

x=np.sin(2*np.pi*5*t)+np.random.normal(0,0.5,t.shape)

#計(jì)算自相關(guān)函數(shù)

acf=sm.tsa.acf(x,nlags=100)

#計(jì)算自相關(guān)函數(shù)的傅里葉變換

frequencies=np.fft.fftfreq(len(acf))

spectrum=np.abs(np.fft.fft(acf))

#繪制自相關(guān)函數(shù)譜

plt.figure(figsize=(10,5))

plt.plot(frequencies,spectrum)

plt.title('自相關(guān)函數(shù)譜')

plt.xlabel('頻率')

plt.ylabel('譜密度')

plt.grid(True)

plt.show()3.2.3解釋在這個例子中，我們首先計(jì)算了時間序列的自相關(guān)函數(shù)。然后，我們使用numpy.fft.fft函數(shù)對自相關(guān)函數(shù)進(jìn)行傅里葉變換，得到自相關(guān)函數(shù)譜。最后，我們繪制了自相關(guān)函數(shù)譜，頻率以Hz為單位，譜密度以線性刻度顯示。通過上述兩種方法，我們可以深入理解時間序列中的周期性特征，這對于預(yù)測和模型構(gòu)建具有重要意義。4譜分析在時間序列中的應(yīng)用4.1季節(jié)性檢測4.1.1原理時間序列的季節(jié)性檢測是通過分析序列中是否存在周期性的模式來實(shí)現(xiàn)的。譜分析，特別是傅里葉變換，可以將時間序列從時域轉(zhuǎn)換到頻域，從而揭示出序列中的周期性成分。在頻域中，周期性模式表現(xiàn)為特定頻率下的峰值，這些峰值對應(yīng)于時間序列中的季節(jié)性周期。4.1.2內(nèi)容傅里葉變換：將時間序列轉(zhuǎn)換為頻譜，以識別周期性模式。功率譜密度：計(jì)算頻譜中各頻率的功率，用于檢測季節(jié)性。周期圖：可視化功率譜密度，直觀展示季節(jié)性周期。4.1.3示例代碼假設(shè)我們有一個時間序列數(shù)據(jù)，表示一年中每個月的氣溫變化，我們想要檢測其中的季節(jié)性。importnumpyasnp

importmatplotlib.pyplotasplt

fromscipy.signalimportperiodogram

#創(chuàng)建模擬數(shù)據(jù)

months=np.arange(1,13)

temperatures=np.sin(2*np.pi*months/12)+np.random.normal(0,0.5,12)

#計(jì)算功率譜密度

frequencies,psd=periodogram(temperatures,fs=12)

#繪制周期圖

plt.figure(figsize=(10,5))

plt.semilogy(frequencies,psd)

plt.title('周期圖')

plt.xlabel('頻率')

plt.ylabel('功率譜密度')

plt.grid(True)

plt.show()4.1.4解釋上述代碼中，我們首先創(chuàng)建了一個模擬的氣溫時間序列，其中包含了一個明顯的年度周期性。然后，我們使用scipy.signal.periodogram函數(shù)計(jì)算了該序列的功率譜密度。最后，我們繪制了周期圖，可以看到在頻率為1（對應(yīng)一年周期）處有一個明顯的峰值，這表明數(shù)據(jù)中存在明顯的年度季節(jié)性。4.2周期性識別4.2.1原理周期性識別是通過分析時間序列的頻譜來確定其周期性的過程。這通常涉及到識別功率譜密度圖中的顯著峰值，這些峰值對應(yīng)于時間序列中的周期。周期性識別對于預(yù)測和模型構(gòu)建至關(guān)重要，因?yàn)樗梢詭椭覀兝斫鈹?shù)據(jù)的內(nèi)在模式。4.2.2內(nèi)容頻譜分析：通過傅里葉變換得到頻譜。峰值檢測：在頻譜中識別顯著的峰值。周期計(jì)算：將峰值頻率轉(zhuǎn)換為時間序列的周期。4.2.3示例代碼假設(shè)我們有一個時間序列數(shù)據(jù)，表示某商品的月銷售量，我們想要識別其中的周期性。importnumpyasnp

importmatplotlib.pyplotasplt

fromscipy.signalimportfind_peaks,periodogram

#創(chuàng)建模擬數(shù)據(jù)

months=np.arange(1,13*5+1)

sales=np.sin(2*np.pi*months/12)+np.sin(2*np.pi*months/6)+np.random.normal(0,0.5,len(months))

#計(jì)算功率譜密度

frequencies,psd=periodogram(sales,fs=12)

#繪制周期圖

plt.figure(figsize=(10,5))

plt.semilogy(frequencies,psd)

plt.title('周期圖')

plt.xlabel('頻率')

plt.ylabel('功率譜密度')

plt.grid(True)

#峰值檢測

peaks,_=find_peaks(psd,height=0.1)

#打印周期

forpeakinpeaks:

print(f'檢測到周期：{12/frequencies[peak]:.2f}個月')

plt.plot(frequencies[peaks],psd[peaks],"x")

plt.show()4.2.4解釋在本例中，我們創(chuàng)建了一個包含兩個周期性模式的模擬銷售數(shù)據(jù)：一個周期為12個月，另一個周期為6個月。我們使用scipy.signal.periodogram計(jì)算了銷售數(shù)據(jù)的功率譜密度，并繪制了周期圖。然后，我們使用scipy.signal.find_peaks函數(shù)來檢測頻譜中的顯著峰值。最后，我們計(jì)算了這些峰值對應(yīng)的周期，并在圖上標(biāo)記了這些峰值。從輸出中，我們可以看到檢測到的周期與我們創(chuàng)建數(shù)據(jù)時設(shè)定的周期相匹配，這驗(yàn)證了我們的周期性識別方法的有效性。5譜分析的高級主題5.1譜泄漏與窗函數(shù)5.1.1譜泄漏原理在時間序列的譜分析中，譜泄漏是一個常見問題，尤其在使用快速傅立葉變換(FFT)時。當(dāng)信號的周期與采樣窗口不匹配時，F(xiàn)FT會假設(shè)信號是周期性的，導(dǎo)致信號在頻域的邊緣處產(chǎn)生不連續(xù)性。這種不連續(xù)性在頻域中表現(xiàn)為能量從信號的真實(shí)頻率泄漏到相鄰的頻率，從而降低了譜分辨率和精度。5.1.2窗函數(shù)的使用為了解決譜泄漏問題，通常會在信號上應(yīng)用窗函數(shù)。窗函數(shù)是一種數(shù)學(xué)函數(shù)，用于對信號進(jìn)行加權(quán)，以減少頻譜泄漏。常見的窗函數(shù)包括漢明窗、海明窗、矩形窗、三角窗、高斯窗等。通過在信號的兩端應(yīng)用窗函數(shù)，可以平滑信號的邊緣，減少頻譜泄漏，提高譜分析的準(zhǔn)確性。漢明窗示例importnumpyasnp

importmatplotlib.pyplotasplt

#生成一個正弦波信號

fs=1000#采樣頻率

t=np.arange(0,1,1/fs)#時間向量

f=50#信號頻率

x=np.sin(2*np.pi*f*t)

#應(yīng)用漢明窗

hamming_window=np.hamming(len(x))

x_windowed=x*hamming_window

#計(jì)算FFT

X=np.fft.fft(x)

X_windowed=np.fft.fft(x_windowed)

#計(jì)算頻譜

freq=np.fft.fftfreq(len(x),1/fs)

X_mag=np.abs(X)

X_windowed_mag=np.abs(X_windowed)

#繪制結(jié)果

plt.figure(figsize=(12,6))

plt.plot(freq,X_mag,label='Withoutwindow')

plt.plot(freq,X_windowed_mag,label='WithHammingwindow')

plt.legend()

plt.xlabel('Frequency(Hz)')

plt.ylabel('Magnitude')

plt.title('Spectrumwithandwithoutwindowing')

plt.show()5.1.3窗函數(shù)選擇選擇窗函數(shù)時，需要考慮兩個主要因素：主瓣寬度和旁瓣水平。主瓣寬度決定了譜分辨率，旁瓣水平則影響了譜泄漏的程度。例如，漢明窗和海明窗都有較窄的主瓣和較低的旁瓣，適合于提高譜分辨率和減少泄漏。矩形窗的主瓣最窄，但旁瓣最高，容易產(chǎn)生泄漏。因此，根據(jù)具體的應(yīng)用場景和需求，選擇合適的窗函數(shù)至關(guān)重要。5.2譜分辨率提升技術(shù)5.2.1高分辨率譜分析傳統(tǒng)的FFT方法在頻域的分辨率受到采樣頻率和信號長度的限制。為了提高譜分辨率，可以采用一些高級技術(shù)，如最大熵譜分析(MEM)、多信號分類(MUSIC)、估計(jì)信號參數(shù)的增強(qiáng)型子空間方法(ESPRIT)等。這些方法通常基于信號的自相關(guān)矩陣或協(xié)方差矩陣，通過數(shù)學(xué)模型來估計(jì)信號的頻譜，從而獲得更高的分辨率。最大熵譜分析示例importnumpyasnp

importmatplotlib.pyplotasplt

fromscipy.signalimportperiodogram,spectral

#生成一個正弦波信號

fs=1000#采樣頻率

t=np.arange(0,1,1/fs)#時間向量

f=50#信號頻率

x=np.sin(2*np.pi*f*t)

#使用最大熵譜分析

f,Pxx_den=spectral(x,Fs=fs,window='hamming',sides='onesided',method='maximum_entropy')

#繪制結(jié)果

plt.figure(figsize=(12,6))

plt.semilogy(f,Pxx_den)

plt.xlabel('Frequency(Hz)')

plt.ylabel('Powerspectraldensity')

plt.title('MaximumEntropySpectrum')

plt.show()5.2.2頻譜插值頻譜插值是另一種提高譜分辨率的方法。通過在FFT結(jié)果上應(yīng)用插值算法，可以在原始頻率點(diǎn)之間生成更多的頻譜點(diǎn)，從而提高分辨率。常見的插值方法包括線性插值、三次樣條插值等。然而，頻譜插值并不能真正增加新的信息，它只是在已有的數(shù)據(jù)點(diǎn)之間進(jìn)行估計(jì)，因此在某些情況下可能不會顯著提高分辨率。頻譜插值示例importnumpyasnp

importmatplotlib.pyplotasplt

fromerpolateimportinterp1d

#生成一個正弦波信號

fs=1000#采樣頻率

t=np.arange(0,1,1/fs)#時間向量

f=50#信號頻率

x=np.sin(2*np.pi*f*t)

#計(jì)算FFT

X=np.fft.fft(x)

freq=np.fft.fftfreq(len(x),1/fs)

#頻譜插值

new_freq=np.linspace(freq.min(),freq.max(),1000)

f=interp1d(freq,np.abs(X),kind='cubic')

X_interpolated=f(new_freq)

#繪制結(jié)果

plt.figure(figsize=(12,6))

plt.plot(new_freq,X_interpolated)

plt.xlabel('Frequency(Hz)')

plt.ylabel('Magnitude')

plt.title('InterpolatedSpectrum')

plt.show()5.2.3零填充零填充是另一種提高頻譜分辨率的簡單方法。在信號的末尾添加零值，可以增加信號的長度，從而在頻域中產(chǎn)生更密集的頻率點(diǎn)。然而，零填充同樣不會增加新的信息，它只是通過增加頻率點(diǎn)的密度來提高分辨率的視覺效果。零填充示例importnumpyasnp

importmatplotlib.pyplotasplt

#生成一個正弦波信號

fs=1000#采樣頻率

t=np.arange(0,1,1/fs)#時間向量

f=50#信號頻率

x=np.sin(2*np.pi*f*t)

#零填充

x_padded=np.pad(x,(0,len(x)),mode='constant')

#計(jì)算FFT

X=np.fft.fft(x)

X_padded=np.fft.fft(x_padded)

#計(jì)算頻譜

freq=np.fft.fftfreq(len(x),1/fs)

freq_padded=np.fft.fftfreq(len(x_padded),1/fs)

#繪制結(jié)果

plt.figure(figsize=(12,6))

plt.plot(freq,np.abs(X),label='Original')

plt.plot(freq_padded,np.abs(X_padded),label='Zero-padded')

plt.legend()

plt.xlabel('Frequency(Hz)')

plt.ylabel('Magnitude')

plt.title('Spectrumwithandwithoutzero-padding')

plt.show()通過上述高級主題的探討，我們可以看到，譜泄漏與窗函數(shù)的選擇、譜分辨率的提升技術(shù)如最大熵譜分析、頻譜插值和零填充，都是在時間序列分析中提高譜分析準(zhǔn)確性和分辨率的重要手段。選擇合適的窗函數(shù)和分辨率提升技術(shù)，可以顯著改善時間序列的譜分析結(jié)果。6實(shí)踐案例與工具6.1使用Python進(jìn)行譜分析在時間序列分析中，譜分析是一種用于識別數(shù)據(jù)中周期性模式和頻率成分的統(tǒng)計(jì)方法。Python提供了多種庫，如numpy,scipy,和matplotlib，可以用來執(zhí)行譜分析。下面，我們將通過一個具體的例子來展示如何使用Python進(jìn)行時間序列的譜分析。6.1.1示例：分析正弦波信號假設(shè)我們有一個由兩個不同頻率的正弦波組成的信號，我們想要通過譜分析來識別這些頻率。數(shù)據(jù)生成importnumpyasnp

importmatplotlib.pyplotasplt

fromscipy.signalimportperiodogram

#設(shè)置隨機(jī)種子以確保結(jié)果可復(fù)現(xiàn)

np.random.seed(0)

#生成時間序列數(shù)據(jù)

t=np.linspace(0,10,1000,endpoint=False)

x=np.sin(2*np.pi*5*t)+0.5*np.sin(2*np.pi*20*t)+np.random.normal(0,0.1,t.shape)

#繪制原始信號

plt.figure()

plt.plot(t,x)

plt.title('原始信號')

plt.xlabel('時間')

plt.ylabel('振幅')

plt.show()譜分析使用scipy.signal.periodogram函數(shù)來計(jì)算信號的功率譜密度。#計(jì)算功率譜密度

frequencies,psd=periodogram(x,fs=100)

#繪制功率譜密度

plt.figure()

plt.semilogy(frequencies,psd)

plt.