série temporelle Compute de python PSD

Question

I ont un spectre signaux PSD qui ressemble à:

La gamme de fréquences du PSD est np.linspace (0,2,500). Je veux convertir ce spectre en une série chronologique de 600s. Le code est indiqué ci-dessous:

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 

Cependant, mon résultat final ressemble à:

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

L'axe x est le nombre de référait points de la série chronologique. Cependant, la série chronologique devrait être de 600s. De l'aide? Merci

Answer 1

Le résultat de votre fonction "spectrumToSeries" a la même longueur que le tableau que vous donnez dans le fichier np.fft.ifft. Parce que la fonction ifft retourne un tableau de la même longueur que l'entrée. Donc, parce que votre tableau psdLoad initial a 500 éléments, le tableau "amplitude" a aussi une longueur de 500 éléments, et donc le fichier randomSeries, qui est le résultat de votre fonction.

Je n'ai pas vraiment les différentes entrées de votre fonction. Quel est le premier argument appelé timeSeries? Est-ce une matrice vide de 600 éléments qui attendent le résultat de la fonction?

J'essaie de calculer des séries chronologiques à partir de PSD moi-même, donc j'aimerais voir votre fonction donner un bon résultat!

Je pense que si vous voulez que votre série temporelle soit de 600 éléments, vous devez avoir une "fréquence" et un tableau "psdLoad" de 600 éléments. Donc, ce que j'essaie de faire avec mon ensemble de données est d'adapter mon psdLoad avec une fonction (psdLoad = f (frequency)). Ensuite, je peux définir la taille de mes tableaux à la longueur de la série de temps que je veux à la fin, et de calculer le ifft ...

Mes propres données sont un enregistrement à 1Hz, sur une journée, donc des tableaux de 86400 éléments . Je dois lui appliquer un filtre, en utilisant une méthode avec PSD. Donc, je calcule mon PSD, dont la longueur est de 129 éléments, et une fois que je l'ai filtré, je veux finir avec mes séries temporelles filtrées.

voici mon code:

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

Le tableau Up13_7april, est ma série de temps initial, dans ce code, je suis juste en train de calculer le PSD, puis revenir à une série de temps pour comparer le signal d'origine et le dernier. Voici le résultat:

[Désolé ne peuvent pas poster une image parce que je suis nouveau à stackoverflow]

donc mon processus est de trouver une fonction qui correspond à la PSD. J'utilise la fonction scipy de Python appelée "optimize.curve_fit". Il vous donne juste les meilleurs paramètres pour adapter vos données avec une fonction que vous fournissez. Une fois que j'ai mes paramètres, je crée de nouveaux tableaux PSD et fréquences, de 86400 éléments. Et enfin j'utilise votre fonction "PSDToSeries" pour calculer le timeSeries.

Je suis assez content du résultat ...Je pense que j'ai juste besoin de trouver un meilleur ajustement de mon PSD:

[Désolé ne peut pas poster une image parce que je suis nouveau à stackoverflow]

Toute idée?