从PSD python计算时间序列 [英] Compute time Series from PSD python

问题描述

我有一个看起来像这样的信号频谱PSD:

I have a signal spectrum PSD that looks like :

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

The frequency range of the PSD is np.linspace(0,2,500). I want to convert this spectrum into a time series of 600s . The code is shown below:

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

但是我的最终结果如下:

However my end result looks like:

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

x轴表示时间序列的点数.但是，时间序列应为600s.有什么帮助吗?谢谢

The x axis is refering the number of points of the time series. However, the time series should be of 600s. Any help? Thanks

推荐答案

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

The result of your function "spectrumToSeries" is the same length as the array you give in the np.fft.ifft. Because the ifft function returns an array of the same length as the input. So, because your initial psdLoad array has 500 elements, the "amplitude" array is 500 elements long too, and so as the randomSeries one, which is your function's result.

我并没有真正得到函数的不同输入.第一个称为timeSeries的参数是什么?它是一个有600个元素的空矩阵，等待函数的结果吗?

I don't really get the different inputs of your function. What is the first argument called timeSeries ? Is it an empty matrix of 600 elements awaiting for the result of the function ?

我正在尝试自己从PSD计算时间序列，所以我很希望看到您的函数给出了很好的结果！

I am trying to compute time series from PSD myself so I'd love to see your function give a good result !

我认为，如果您希望时间序列为600个元素，则需要有一个频率"和一个"psdLoad"数组，其中包含600个元素.因此，我要处理的数据集是使psdLoad具有一个函数(psdLoad = f(频率)).然后，我可以将数组的大小设置为最后想要的时间序列的长度，然后计算ifft ...

I think that if you want your time series to be 600 elements, you need to have a "frequency" and a "psdLoad" array of 600 elements. So what I am trying to do with my set of data is to fit my psdLoad with a function (psdLoad = f (frequency)). Then I can set the size of my arrays to the length of the timeseries I want at the end, and compute the ifft...

我自己的数据是一天中以1Hz记录的，因此有86400个元素的数组.我必须使用带有PSD的方法对其应用过滤器.因此，我计算了PSD的长度，即129个元素，对它进行过滤后，我想以过滤后的时间序列结束.

My own data is a record at 1Hz, over a day, so arrays of 86400 elements. I have to apply a filter to it, using a method with PSD. So I compute my PSD, which length is 129 elements, and once I have filtered it I want to end up with my filtered time series.

这是我的代码:

######################################################################"
## 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，然后返回到时间序列以比较原始信号和最终信号.结果如下:

The array Up13_7april, is my initial time series, in this code I am just trying to compute the PSD and then come back to a time serie to compare the original signal and the final one. Here is the result :

[对不起，因为我是stackoverflow的新手，所以无法发布任何图片

[Sorry can't post any picture because I'm new to stackoverflow]

所以我的过程是找到适合PSD的功能.我使用称为"optimize.curve_fit"的Python scipy函数.它只是为您提供最佳参数，以使用您提供的功能来适合您的数据. 一旦有了参数，就可以创建包含86400个元素的新PSD和频率阵列.最后，我使用您的"PSDToSeries"函数来计算timeSeries.

So my process is to find a function that fits the PSD. I use the Python scipy function called "optimize.curve_fit". It just gives you the best parameters to fit your data with a function that you provide. Once I have my parameters, I create new PSD and frequency arrays, of 86400 elements. And finally I use your "PSDToSeries" function to compute the timeSeries.

我对结果感到非常满意……我想我只需要找到更适合我的PSD即可:

I'm quite happy with the result... I think I just need to find a better fit of my PSD :

[对不起，因为我是stackoverflow的新手，所以无法发布任何图片

[Sorry can't post any picture because I'm new to stackoverflow]

有什么主意吗?

这篇关于从PSD python计算时间序列的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！