从 PSD python 计算时间序列答案

【问题标题】：Compute time Series from PSD python从 PSD python 计算时间序列
【发布时间】：2015-03-17 09:23:13
【问题描述】：

我的信号频谱 PSD 看起来像：

PSD 的频率范围是 np.linspace(0,2,500)。我想将此频谱转换为 600s 的时间序列。代码如下：

def spectrumToSeries(timeSeries,frequency,psdLoad):
    ''' 
    Function that gicen a PSD converts into a time series

    '''
    #
    #Obtian interval frequency
    df=frequency[2]-frequency[1]    

    #Obtian the spectrum amplitudes
    amplitude=np.sqrt(2*np.array(psdLoad)*df)

    #Pre allocation of matrices
    epsilon=np.zeros((len(amplitude)))
    randomSeries=np.zeros((len(amplitude)))


    #Create time series from spectrum
    #Generate random phases between [-2pi,2pi]
    epsilon=-np.pi + 2*np.pi*np.random.randn(1,len(amplitude))

    #Inverse Fourier
    randomSeries=len(timeSeries)*np.real(np.fft.ifft(amplitude*np.exp(epsilon*1j*2*np.pi))));

    return randomSeries

但是我的最终结果是这样的：

timeSeries = spectrumToSeries(thrustBladed,param.frequency,analyticalThrustPSD[iwind])

x 轴表示时间序列的点数。但是，时间序列应该是 600 秒。有什么帮助吗？谢谢

【问题讨论】：

标签： python time-series fft spectrum

【解决方案1】：

您的函数“spectrumToSeries”的结果与您在 np.fft.ifft 中给出的数组长度相同。因为 ifft 函数返回一个与输入长度相同的数组。因此，因为您的初始 psdLoad 数组有 500 个元素，所以“幅度”数组也是 500 个元素长，randomSeries 也是如此，这是您的函数的结果。

我并没有真正了解您函数的不同输入。第一个称为 timeSeries 的参数是什么？它是等待函数结果的 600 个元素的空矩阵吗？

我正在尝试自己从 PSD 计算时间序列，所以我很想看到你的函数给出一个好的结果！

我认为，如果您希望您的时间序列为 600 个元素，您需要有一个包含 600 个元素的“频率”和一个“psdLoad”数组。所以我试图用我的数据集做的是用一个函数（psdLoad = f（频率））来拟合我的 psdLoad。然后我可以将数组的大小设置为最后我想要的时间序列的长度，并计算 ifft...

我自己的数据是一天内 1Hz 的记录，因此是 86400 个元素的数组。我必须使用带有 PSD 的方法对其应用过滤器。所以我计算了我的 PSD，它的长度是 129 个元素，一旦我过滤了它，我想最终得到我过滤的时间序列。

这是我的代码：

######################################################################"
## Computation of spectrum values : PSD & frequency ##
######################################################################"

psd_ampl0, freq = mlab.psd(Up13_7april, NFFT=256, Fs=1, detrend=mlab.detrend_linear, window=mlab.window_hanning, noverlap=0.5, sides='onesided')

################################################"
## Computation of the time series from the PSD ##
################################################"


def PSDToSeries(lenTimeSeries,freq,psdLoad):
    ''' 
    Function that gicen a PSD converts into a time series

    '''
    #
    #Obtian interval frequency
    df=freq[2]-freq[1]    
    print('df = ', df)

    #Obtian the spectrum amplitudes
    amplitude=(2*psdLoad*df)**0.5

    #Pre allocation of matrices
    epsilon=np.zeros((len(amplitude)))
    randomSeries=np.zeros((len(amplitude)))


    #Create time series from spectrum
    #Generate random phases between [-2pi,2pi]
    epsilon=-np.pi + 2*np.pi*np.random.randn(1,len(amplitude))

    #Inverse Fourier
    randomSeries=lenTimeSeries*np.real(np.fft.ifft(amplitude*np.exp(epsilon*1j*2*np.pi)));

    return randomSeries

#-------------------------------------------------------------------------

#########################################################"
## Fitting a function on the PSD to add it more points ##
#########################################################"

#def fitting_function(freq,a,b,c,d,e,f):
    #return a*(freq**5)+b*(freq**4)+c*(freq**3)+d*(freq**2)+e*freq+f

def fitting_function(freq,a,b,c):
    return a*np.exp(freq*b)

# Look for the best fitting parameters of the choosen fitting function #

param_opt, pcov = optim.curve_fit(fitting_function,freq[1:],psd_ampl0[1:])

print('The best fitting parameters are : ',param_opt)

# Definition of the new PSD and frequency arrays extended to 86400 elements #

freq_extend = np.linspace(min(freq),max(freq), 86400)

psd_extend = fitting_function(freq_extend,param_opt[0], param_opt[1], param_opt[2])

#print(psd_allonge)

ts_length = Up13_7april.shape[0] #Length of the timeSeries I want to compute

print('ts_length = ', ts_length)

tsFromPSD = PSDToSeries(ts_length, freq_allonge, psd_allonge)

print('shape tsFromPSD : ', tsFromPSD.shape)


##################"
## Plot section ##
##################"

plt.figure(1)
plt.plot(freq[1:] ,psd_ampl0[1:],marker=',', ls='-',color='SkyBlue', label='original PSD')
plt.plot(freq_allonge, psd_allonge,  marker=',', ls='-',color='DarkGreen', label='PSD rallonge')
plt.xlabel('Frequency [Hz]')
plt.ylabel('PSD of raw velocity module [(m/s)²/Hz]')
plt.grid(True)
plt.legend()


plt.figure(2)
plt.plot_date(time7april,Up13_7april, xdate=True, ydate=False, marker=',', ls='-', c='Grey', label='Original Signal')
plt.plot_date(time7april, tsFromPSD[0],xdate=True, ydate=False, marker=',', ls='-', label='After inverse PSD')
plt.suptitle('Original and Corrected time series for the 7th of April')
plt.grid(True)
plt.legend()

plt.show()

数组 Up13_7april 是我的初始时间序列，在这段代码中，我只是尝试计算 PSD，然后返回时间序列来比较原始信号和最终信号。结果如下：

[抱歉无法发布任何图片，因为我是stackoverflow的新手]

所以我的过程是找到一个适合 PSD 的函数。我使用名为“optimize.curve_fit”的 Python scipy 函数。它只是为您提供最佳参数，以使您的数据与您提供的功能相匹配。获得参数后，我将创建 86400 个元素的新 PSD 和频率阵列。最后，我使用您的“PSDToSeries”函数来计算 timeSeries。

我对结果很满意...我想我只需要找到更适合我的 PSD ：

[抱歉无法发布任何图片，因为我是stackoverflow的新手]

有什么想法吗？

【讨论】：

是的，时间序列是一个空向量！那么我的脚本中的错误可能是什么？
我认为使用 np.fft.ifft 会发生什么，您计算一个大小为 500 的数组。然后，您只需将其乘以 600，因为我假设“timeSeries”大小为 600。但最后你没有得到一个包含 600 个元素的数组。我认为如果你希望你的时间序列是 600 个元素，你需要有一个 600 个元素的“频率”和一个“psdLoad”数组。所以我试图用我的数据集做的是用一个函数（psdLoad = f（频率））来拟合我的 psdLoad。然后我可以将数组的大小设置为最后我想要的时间序列的长度，并计算 ifft...