J'ai du mal à générer un ensemble de séries temporelles colorées stationnaires, étant donné leur matrice de covariance (leurs densités spectrales de puissance (PSD) et leurs densités spectrales de puissance croisée (CSD)).
Je sais que, compte tenu de deux séries chronologiques et , je peux estimer les densités spectrales de puissance (STDCP) et des densités spectrales croisées (SDR) en utilisant de nombreuses routines largement disponibles, comme psd()
et csd()
fonctions Matlab , etc. Les PSD et les CSD constituent la matrice de covariance:
Que se passe-t-il si je veux faire l'inverse? Compte tenu de la matrice de covariance, comment générer une réalisation de et ?
Veuillez inclure toute théorie de base ou signaler tous les outils existants qui font cela (tout ce qui en Python serait formidable).
Ma tentative
Vous trouverez ci-dessous une description de ce que j'ai essayé et des problèmes que j'ai constatés. C'est un peu long à lire, et désolé s'il contient des termes qui ont été mal utilisés. Si ce qui est erroné peut être signalé, ce serait très utile. Mais ma question est celle en gras ci-dessus.
- Les PSD et les CSD peuvent être écrits comme la valeur attendue (ou moyenne d'ensemble) des produits des transformées de Fourier de la série chronologique. Ainsi, la matrice de covariance peut s'écrire:
où
- Une matrice de covariance est une matrice hermitienne, ayant des valeurs propres réelles qui sont soit nulles soit positives. Donc, il peut être décomposé en
oùest une matrice diagonale dont les éléments non nuls sont les racines carréesdes valeurs propresde; est la matrice dont les colonnes sont les vecteurs propres orthonormés de; est la matrice d'identité.
- La matrice d'identité est écrite comme
oùet sont des séries de fréquences non corrélées et complexes avec une moyenne et une variance unitaire nulles.
- En utilisant 3. dans 2., puis comparer avec 1. Les transformées de Fourier de la série chronologique sont:
- Les séries chronologiques peuvent ensuite être obtenues en utilisant des routines comme la transformée de Fourier rapide inverse.
J'ai écrit une routine en Python pour ce faire:
def get_noise_freq_domain_CovarMatrix( comatrix , df , inittime , parityN , seed='none' , N_previous_draws=0 ) :
"""
returns the noise time-series given their covariance matrix
INPUT:
comatrix --- covariance matrix, Nts x Nts x Nf numpy array
( Nts = number of time-series. Nf number of positive and non-Nyquist frequencies )
df --- frequency resolution
inittime --- initial time of the noise time-series
parityN --- is the length of the time-series 'Odd' or 'Even'
seed --- seed for the random number generator
N_previous_draws --- number of random number draws to discard first
OUPUT:
t --- time [s]
n --- noise time-series, Nts x N numpy array
"""
if len( comatrix.shape ) != 3 :
raise InputError , 'Input Covariance matrices must be a 3-D numpy array!'
if comatrix.shape[0] != comatrix.shape[1] :
raise InputError , 'Covariance matrix must be square at each frequency!'
Nts , Nf = comatrix.shape[0] , comatrix.shape[2]
if parityN == 'Odd' :
N = 2 * Nf + 1
elif parityN == 'Even' :
N = 2 * ( Nf + 1 )
else :
raise InputError , "parityN must be either 'Odd' or 'Even'!"
stime = 1 / ( N*df )
t = inittime + stime * np.arange( N )
if seed == 'none' :
print 'Not setting the seed for np.random.standard_normal()'
pass
elif seed == 'random' :
np.random.seed( None )
else :
np.random.seed( int( seed ) )
print N_previous_draws
np.random.standard_normal( N_previous_draws ) ;
zs = np.array( [ ( np.random.standard_normal((Nf,)) + 1j * np.random.standard_normal((Nf,)) ) / np.sqrt(2)
for i in range( Nts ) ] )
ntilde_p = np.zeros( ( Nts , Nf ) , dtype=complex )
for k in range( Nf ) :
C = comatrix[ :,:,k ]
if not np.allclose( C , np.conj( np.transpose( C ) ) ) :
print "Covariance matrix NOT Hermitian! Unphysical."
w , V = sp_linalg.eigh( C )
for m in range( w.shape[0] ) :
w[m] = np.real( w[m] )
if np.abs(w[m]) / np.max(w) < 1e-10 :
w[m] = 0
if w[m] < 0 :
print 'Negative eigenvalue! Simulating unpysical signal...'
ntilde_p[ :,k ] = np.conj( np.sqrt( N / (2*stime) ) * np.dot( V , np.dot( np.sqrt( np.diag( w ) ) , zs[ :,k ] ) ) )
zerofill = np.zeros( ( Nts , 1 ) )
if N % 2 == 0 :
ntilde = np.concatenate( ( zerofill , ntilde_p , zerofill , np.conj(np.fliplr(ntilde_p)) ) , axis = 1 )
else :
ntilde = np.concatenate( ( zerofill , ntilde_p , np.conj(np.fliplr(ntilde_p)) ) , axis = 1 )
n = np.real( sp.ifft( ntilde , axis = 1 ) )
return t , n
J'ai appliqué cette routine aux PSD et aux CSD, dont les expressions analytiques ont été obtenues à partir de la modélisation d'un détecteur avec lequel je travaille. L'important est qu'à toutes les fréquences, ils constituent une matrice de covariance (enfin au moins ils passent toutes ces if
déclarations dans la routine). La matrice de covariance est 3x3. Les 3 séries chronologiques ont été générées environ 9000 fois, et les PSD et CSD estimés, moyennés sur toutes ces réalisations sont tracés ci-dessous avec les analyses. Bien que les formes globales concordent, il y a des caractéristiques bruyantes notables à certaines fréquences dans les CSD (Fig.2). Après un gros plan sur les pics des PSD (Fig.3), j'ai remarqué que les PSD sont en fait sous-estimés, et que les caractéristiques bruyantes dans les CSD se produisent à peu près aux mêmes fréquences que les pics dans les PSD. Je ne pense pas que ce soit une coïncidence, et que d'une manière ou d'une autre, une fuite de puissance des PSD vers les CSD. Je me serais attendu à ce que les courbes se superposent, avec autant de réalisations des données.
la source
Réponses:
Étant donné que vos signaux sont stationnaires, une approche simple serait d'utiliser le bruit blanc comme base et de le filtrer pour l'adapter à vos PSD. Une façon de calculer ces coefficients de filtre consiste à utiliser la prédiction linéaire .
Il semble qu'il y ait une fonction python pour cela, essayez-la:
Si vous le souhaitez (je n'ai utilisé que l'équivalent MATLAB). Il s'agit d'une approche utilisée dans le traitement de la parole, où les formants sont estimés de cette façon.
la source
Un peu tard pour la fête, comme d'habitude, mais je vois une activité récente donc je vais mes deux yens.
Tout d'abord, je ne peux pas blâmer la tentative de PO - cela me semble juste. Les écarts pourraient être dus à des problèmes avec des échantillons finis, par exemple un biais positif de l'estimation de la puissance du signal.
Cependant, je pense qu'il existe des moyens plus simples de générer des séries temporelles à partir de la matrice de densité spectrale croisée (CPSD, c'est ce que l'OP a appelé matrice de covariance).
Une approche paramétrique consiste à utiliser le CPSD pour obtenir une description autorégressive, puis à l'utiliser pour générer la série chronologique. Dans matlab, vous pouvez le faire en utilisant les outils de causalité de Granger (par exemple la boîte à outils de causalité Multivaraite Granger, Seth, Barnett ). La boîte à outils est très simple à utiliser. Étant donné que l'existence du CPSD garantit une description autorégressive, cette approche est exacte. (pour plus d'informations sur le CPSD et l'autorégression, voir "Mesure de la dépendance linéaire et du feedback entre plusieurs séries chronologiques" par Geweke, 1982, ou de nombreux articles d'Aneil Seth + Lionel Barnett, pour obtenir l'image complète).
Potentiellement plus simple, notant que le CPSD peut être formé en appliquant le fft à la covariance automatique (donnant la diagonale du CPSD, c'est-à-dire la puissance des signaux) et la covariance croisée (donnant les éléments diagonaux éteints, c'est-à-dire la puissance croisée). Ainsi, en appliquant l'inverse fft au CPSD, nous pouvons obtenir l'autocorrélation et l'auto covariance. Nous pouvons ensuite les utiliser pour générer des échantillons de nos données.
J'espère que cela t'aides. Veuillez laisser toutes les demandes d'informations dans les commentaires et je vais essayer de répondre.
la source