Simulación de series temporales de potencia dada y densidades espectrales cruzadas

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 , puedo estimar sus densidades espectrales de potencia (PSD) y densidades espectrales cruzadas (CSD) utilizando muchas rutinas ampliamente disponibles, como y funciones en Matlab, etc. Los PSD y CSD forman la matriz de covarianza: $y_{I}(t)$$y_{J}(t)$psd()csd()

$$\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 una función de frecuencia . $f$

¿Qué pasa si quiero hacer lo contrario? Dada la matriz de covarianza, ¿cómo genero una realización de y ?$y_{I}(t)$$y_{J}(t)$

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 se muestra una descripción de lo que he intentado y los problemas que he notado. Es una lectura un poco larga, y lo siento si contiene términos que han sido mal utilizados. Si se puede señalar lo que es erróneo, sería muy útil. Pero mi pregunta es la que está en negrita arriba.

1. Los PSD y CSD pueden escribirse como el valor esperado (o promedio de conjunto) de los productos de las transformadas de Fourier de la serie temporal. Entonces, la matriz de covarianza se puede escribir como: donde $$\mathbf{C}(f) = \frac{2}{\tau} \langle \mathbf{Y}^{\dagger}(f) \mathbf{Y}(f) \rangle \;,$$ $$\mathbf{Y}(f) = \left( \begin{array}{cc} \tilde{y}_{I}(f) & \tilde{y}_{J}(f) \end{array} \right) \;.$$
2. Una matriz de covarianza es una matriz hermitiana, que tiene valores propios reales que son cero o positivos. Por lo tanto, se puede descomponer en donde es una matriz diagonal cuyos elementos distintos de cero son las raíces cuadradas de los valores propios de ; es la matriz cuyas columnas son los vectores propios ortonormales de ; $$\mathbf{C}(f) = \mathbf{X}(f) \boldsymbol\lambda^{\frac{1}{2}}(f) \: \mathbf{I} \: \boldsymbol\lambda^{\frac{1}{2}}(f) \mathbf{X}^{\dagger}(f) \;,$$ $\lambda^{\frac{1}{2}}(f)$$\mathbf{C}(f)$$\mathbf{X}(f)$$\mathbf{C}(f)$$\mathbf{I}$ es la matriz de identidad.
3. La matriz de identidad se escribe como donde y son series de frecuencia complejas y sin correlación con media cero y varianza unitaria. $$\mathbf{I} = \langle \mathbf{z}^{\dagger}(f) \mathbf{z}(f) \rangle \;,$$ $$\mathbf{z}(f) = \left( \begin{array}{cc} z_{I}(f) & z_{J}(f) \end{array} \right) \;,$$ $\{z_{i}(f)\}_{i=I,J}$
4. Al usar 3. en 2., y luego comparar con 1. Las transformadas de Fourier de la serie temporal son: $$\mathbf{Y}(f) = \sqrt{ \frac{\tau}{2} } \mathbf{z}(f) \: \boldsymbol\lambda^{\frac{1}{2}}(f) \: \mathbf{X}^{\dagger}(f) \;.$$
5. La serie temporal se puede obtener utilizando rutinas como la transformada rápida inversa de Fourier.

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 , n

He aplicado esta rutina a PSD y CSD, cuyas expresiones analíticas se han obtenido del modelado de algún detector con el que estoy trabajando. Lo importante es que en todas las frecuencias, forman una matriz de covarianza (bueno, al menos pasan todas esas ifdeclaraciones en la rutina). La matriz de covarianza es 3x3. Las 3 series temporales se han generado aproximadamente 9000 veces, y los PSD y CSD estimados, promediados sobre todas estas realizaciones, se representan a continuación con los analíticos. Si bien las formas generales están de acuerdo, hay características ruidosas notables en ciertas frecuencias en los CSD (Fig.2). Después de un primer plano alrededor de los picos en los PSD (Fig.3), noté que los PSD en realidad están subestimados, y que las características ruidosas en los CSD ocurren casi a las mismas frecuencias que los picos en los PSD. No creo que esto sea una coincidencia, y que de alguna manera el poder se está filtrando de los PSD a los CSD. Hubiera esperado que las curvas estuvieran una encima de la otra, con tantas realizaciones de los datos.

Dado que sus señales son estacionarias, un enfoque simple sería utilizar el ruido blanco como base y filtrarlo para que se ajuste a sus PSD. Una forma de calcular estos coeficientes de filtro es usar predicción lineal .

Parece que hay una función de Python para ello, pruébalo:

from scikits.talkbox import lpc

Si lo desea (solo he usado el equivalente de MATLAB). Este es un enfoque utilizado en el procesamiento del habla, donde los formantes se estiman de esta manera.

Un poco tarde para la fiesta, como de costumbre, pero veo algo de actividad reciente, así que voy a mis dos yenes.

Primero, no puedo criticar el intento de los OP: me parece correcto. Las discrepancias podrían deberse a problemas con las muestras finitas, por ejemplo, el sesgo positivo de la estimación de la potencia de la señal.

Sin embargo, creo que hay formas más simples de generar series de tiempo a partir de la matriz de densidad espectral cruzada (CPSD, esto es lo que el OP llamó matriz de covarianza).

Un enfoque paramétrico es usar el CPSD para obtener una descripción autorregresiva y luego usarla para generar la serie temporal. En matlab puede hacerlo utilizando las herramientas de causalidad de Granger (por ejemplo , la caja de herramientas de causalidad Graniva de Multivaraite, Seth, Barnett ). La caja de herramientas es muy fácil de usar. Dado que la existencia del CPSD garantiza una descripción autorregresiva, este enfoque es exacto. (para obtener más información sobre el CPSD y la autorregresión, consulte "Medición de la dependencia lineal y la retroalimentación entre series temporales múltiples" de Geweke, 1982, o muchos de los documentos de Aneil Seth + Lionel Barnett, para obtener una imagen completa).

Potencialmente más simple, señala que el CPSD se puede formar aplicando el fft a la covarianza automática (que da la diagonal de CPSD, es decir, la potencia de las señales) y la covarianza cruzada (que da los elementos diagonales, es decir, la potencia cruzada). Por lo tanto, aplicando el fft inverso al CPSD podemos obtener la autocorrelación y la covarianza automática. Luego podemos usarlos para generar muestras de nuestros datos.

Espero que esto ayude. Deje cualquier solicitud de información en los comentarios y trataré de responder.