Simulation von Zeitreihen bei gegebener Leistung und Kreuzspektraldichte

Aufgrund ihrer Kovarianzmatrix (ihre Leistungsspektraldichten (PSDs) und Kreuzleistungsspektraldichten (CSDs)) habe ich Probleme, eine Reihe stationärer farbiger Zeitreihen zu erstellen.

Ich weiß, dass ich mit zwei Zeitreihen und ihre Leistungsspektraldichten (PSDs) und Kreuzspektraldichten (CSDs) unter Verwendung vieler allgemein verfügbarer Routinen abschätzen kann, wie z und Funktionen in Matlab usw. Die PSDs und CSDs bilden die Kovarianzmatrix: $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)\;,$$ was im Allgemeinen eine Funktion der Frequenz . $f$

Was passiert, wenn ich das Gegenteil tun möchte? Wie kann ich unter Berücksichtigung der Kovarianzmatrix eine Realisierung von und erzeugen ?$y_{I}(t)$$y_{J}(t)$

Bitte fügen Sie eine Hintergrundtheorie hinzu oder weisen Sie auf vorhandene Tools hin, die dies tun (alles in Python wäre großartig).

Mein Versuch

Nachstehend finden Sie eine Beschreibung meiner Versuche und der Probleme, die mir aufgefallen sind. Es ist ein bisschen langwierig und es tut mir leid, wenn es Begriffe enthält, die missbraucht wurden. Wenn darauf hingewiesen werden kann, was fehlerhaft ist, wäre dies sehr hilfreich. Aber meine Frage ist die oben fett gedruckte.

1. Die PSDs und CSDs können als Erwartungswert (oder Ensemble-Durchschnitt) der Produkte der Fourier-Transformationen der Zeitreihen geschrieben werden. Die Kovarianzmatrix kann also wie geschrieben werden: wo $$\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. Eine Kovarianzmatrix ist eine Hermitianische Matrix mit reellen Eigenwerten, die entweder Null oder positiv sind. Es kann also zerlegt werden in wobei ist eine Diagonalmatrix, deren Nicht-Null-Elemente die Quadratwurzeln der Eigenwerte von ; ist die Matrix, deren Spalten die orthonormalen Eigenvektoren von ; $$\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}$ ist die Identitätsmatrix.
3. Die Identitätsmatrix wird geschrieben als wobei und sind unkorrelierte und komplexe Frequenzreihen mit einem Mittelwert von Null und einer Einheitsvarianz. $$\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. Durch Verwenden von 3. in 2. und dann Vergleichen mit 1. Die Fourier-Transformationen der Zeitreihe lauten: $$\mathbf{Y}(f) = \sqrt{ \frac{\tau}{2} } \mathbf{z}(f) \: \boldsymbol\lambda^{\frac{1}{2}}(f) \: \mathbf{X}^{\dagger}(f) \;.$$
5. Die Zeitreihen können dann unter Verwendung von Routinen wie der inversen schnellen Fourier-Transformation erhalten werden.

Ich habe dazu in Python eine Routine geschrieben:

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

Ich habe diese Routine auf PSDs und CSDs angewendet, deren analytische Ausdrücke aus der Modellierung eines Detektors stammen, mit dem ich arbeite. Wichtig ist, dass sie bei allen Frequenzen eine Kovarianzmatrix bilden (zumindest bestehen sie alle diese ifAnweisungen in der Routine). Die Kovarianzmatrix ist 3x3. Die 3 Zeitreihen wurden ungefähr 9000-mal generiert, und die geschätzten PSDs und CSDs, gemittelt über alle diese Realisierungen, sind unten mit den analytischen dargestellt. Während die Gesamtformen einverstanden ist , gibt es spürbare laut Funktionen bei bestimmten Frequenzen in den Zentralverwahrern (Abb.2). Nach einer Nahaufnahme der Peaks in den PSDs (Abb. 3) bemerkte ich, dass die PSDs tatsächlich unterschätzt werdenund dass die Rauschmerkmale in den CSDs bei ungefähr den gleichen Frequenzen auftreten wie die Spitzen in den PSDs. Ich denke nicht, dass dies ein Zufall ist und dass irgendwie Strom aus den PSDs in die CSDs fließt. Ich hätte erwartet, dass die Kurven bei so vielen Erkenntnissen der Daten übereinander liegen.

Da Ihre Signale stationär sind, besteht ein einfacher Ansatz darin, weißes Rauschen als Basis zu verwenden und es zu filtern, um es an Ihre PSDs anzupassen. Eine Möglichkeit zur Berechnung dieser Filterkoeffizienten ist die Verwendung einer linearen Vorhersage .

Es scheint, dass es eine Python-Funktion dafür gibt. Probieren Sie es aus:

from scikits.talkbox import lpc

Wenn Sie möchten (ich habe nur das MATLAB-Äquivalent verwendet). Dies ist ein in der Sprachverarbeitung verwendeter Ansatz, bei dem Formanten auf diese Weise geschätzt werden.

Wie immer ein bisschen zu spät zur Party, aber ich sehe einige neue Aktivitäten, also werde ich meine zwei Yen.

Erstens kann ich den OP-Versuch nicht beanstanden - er sieht für mich richtig aus. Die Abweichungen könnten auf Probleme mit endlichen Abtastwerten zurückzuführen sein, beispielsweise auf eine positive Vorspannung der Signalleistungsschätzung.

Ich denke jedoch, dass es einfachere Möglichkeiten gibt, Zeitreihen aus der Kreuzspektraldichtematrix (CPSD, das ist, was das OP als Kovarianzmatrix bezeichnet) zu erzeugen.

Ein parametrischer Ansatz besteht darin, die CPSD zu verwenden, um eine autoregressive Beschreibung zu erhalten und diese dann zum Generieren der Zeitreihen zu verwenden. In matlab können Sie dies mit den Granger-Werkzeugen für die Kausalität (z. B. Multivaraite Granger Causality Toolbox, Seth, Barnett ) tun . Die Toolbox ist sehr einfach zu bedienen. Da die Existenz der CPSD eine autoregressive Beschreibung garantiert, ist dieser Ansatz genau. (Weitere Informationen zu CPSD und Autoregression finden Sie in "Messung der linearen Abhängigkeit und Rückkopplung zwischen mehreren Zeitreihen" von Geweke, 1982, oder in vielen Artikeln von Aneil Seth + Lionel Barnett, um ein vollständiges Bild zu erhalten.)

Potenziell einfacher ist die Feststellung, dass die CPSD gebildet werden kann, indem fft auf die Autokovarianz (Angabe der Diagonale von CPSD, dh der Signalleistung) und die Kreuzkovarianz (Angabe der Off-Diagonal-Elemente, dh der Kreuzleistung) angewendet wird. Somit können wir durch Anwenden des inversen fft auf die CPSD die Autokorrelation und Autokovarianz erhalten. Wir können diese dann verwenden, um Proben unserer Daten zu generieren.

Hoffe das hilft. Bitte hinterlassen Sie alle Anfragen für Informationen in den Kommentaren und ich werde versuchen, zu adressieren.