Tengo problemas para generar un conjunto de series temporales de colores estacionarios, dada su matriz de covarianza (sus densidades espectrales de potencia (PSD) y sus densidades espectrales de potencia cruzada (CSD)).
Sé que, dadas dos series temporales $y_{I}(t)$ et $y_{J}(t)$ puedo estimar sus densidades espectrales de potencia (PSD) y sus densidades espectrales cruzadas (CSD) utilizando muchas rutinas ampliamente disponibles, tales como psd() et csd() funciones en Matlab, etc. Los PSD y CSD forman la matriz de covarianza: $$ \mathbf{C}(f) = \left( \begin{array}{cc} P_{II}(f) & P_{IJ}(f)\\ P_{JI}(f) & P_{JJ}(f) \end{array} \right)\;, $$ que en general es función de la frecuencia $f$ .
¿Qué pasa si quiero hacer lo contrario? Dada la matriz de covarianzas, ¿cómo puedo generar una realización de $y_{I}(t)$ et $y_{J}(t)$ ?
Por favor, incluya cualquier teoría de fondo, o señale cualquier herramienta existente que haga esto (cualquier cosa en Python sería genial).
Mi intento
A continuación describo lo que he probado y los problemas que he observado. Es un poco largo de leer, y lo siento si contiene términos que han sido mal utilizados. Si se puede señalar lo que es erróneo, sería de gran ayuda. Pero mi pregunta es la que está en negrita arriba.
- Las PSD y CSD se pueden escribir como el valor de la expectativa (o media del conjunto) de los productos de las transformadas de Fourier de las series temporales. Así, la matriz de covarianza puede escribirse como: $$ \mathbf{C}(f) = \frac{2}{\tau} \langle \mathbf{Y}^{\dagger}(f) \mathbf{Y}(f) \rangle \;, $$ donde $$ \mathbf{Y}(f) = \left( \begin{array}{cc} \tilde{y}_{I}(f) & \tilde{y}_{J}(f) \end{array} \right) \;. $$
- Una matriz de covarianza es una matriz hermitiana, con valores propios reales que son cero o positivos. Por lo tanto, se puede descomponer en $$ \mathbf{C}(f) = \mathbf{X}(f) \boldsymbol\lambda^{\frac{1}{2}}(f) \: \mathbf{I} \: \boldsymbol\lambda^{\frac{1}{2}}(f) \mathbf{X}^{\dagger}(f) \;, $$ donde $\lambda^{\frac{1}{2}}(f)$ es una matriz diagonal cuyos elementos distintos de cero son las raíces cuadradas de $\mathbf{C}(f)$ valores propios; $\mathbf{X}(f)$ es la matriz cuyas columnas son los vectores propios ortonormales de $\mathbf{C}(f)$ ; $\mathbf{I}$ es la matriz de identidad.
- La matriz identidad se escribe como $$ \mathbf{I} = \langle \mathbf{z}^{\dagger}(f) \mathbf{z}(f) \rangle \;, $$ donde $$ \mathbf{z}(f) = \left( \begin{array}{cc} z_{I}(f) & z_{J}(f) \end{array} \right) \;, $$ et $\{z_{i}(f)\}_{i=I,J}$ son series de frecuencias no correlacionadas y complejas con media cero y varianza unitaria.
- Utilizando 3. en 2., y comparando después con 1. Las transformadas de Fourier de las series temporales son: $$ \mathbf{Y}(f) = \sqrt{ \frac{\tau}{2} } \mathbf{z}(f) \: \boldsymbol\lambda^{\frac{1}{2}}(f) \: \mathbf{X}^{\dagger}(f) \;. $$
- A continuación, las series temporales pueden obtenerse mediante rutinas como la transformada rápida de Fourier inversa.
He escrito una rutina en Python para hacer esto:
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 , nHe aplicado esta rutina a PSDs y CSDs, cuyas expresiones analíticas se han obtenido a partir del modelado de algún detector con el que estoy trabajando. Lo importante es que en todas las frecuencias conforman una matriz de covarianza (bueno al menos pasan todas esas if en la rutina). La matriz de covarianza es 3x3. Las 3 series temporales se han generado unas 9.000 veces, y las PSD y CSD estimadas, promediadas sobre todas estas realizaciones, se representan a continuación con las analíticas. Aunque las formas generales coinciden, se observan características ruidosas en determinadas frecuencias de los CSD (Fig.2). Tras un acercamiento a los picos de las PSD (Fig.3), observé que las PSD son en realidad subestimado y que las características ruidosas de los CSD se producen casi a las mismas frecuencias que los picos de los PSD. No creo que sea una coincidencia y que, de alguna manera, se esté filtrando energía de los PSD a los CSD. Habría esperado que las curvas se superpusieran con tantas realizaciones de los datos.
