Symulowanie szeregów czasowych przy danej mocy i gęstościach widm poprzecznych

Mam problem z wygenerowaniem zestawu stacjonarnych kolorowych szeregów czasowych, biorąc pod uwagę ich macierz kowariancji (ich gęstości widmowe mocy (PSD) i gęstości widmowe mocy krzyżowej (CSD)).

Wiem, że biorąc pod uwagę dwie serie czasowe i , mogę oszacować ich gęstość widmową mocy (PSD) i gęstość krzyżową widmową (CSD) przy użyciu wielu powszechnie dostępnych procedur, takich jak i funkcje w Matlabie itp. PSD i CSD tworzą macierz kowariancji: $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)\;,$$ f który jest generalnie funkcją częstotliwości . $f$

Co się stanie, jeśli chcę zrobić odwrotnie? Biorąc pod uwagę macierz kowariancji, jak wygenerować realizację i ?$y_{I}(t)$$y_{J}(t)$

Podaj dowolną teorię tła lub wskaż istniejące narzędzia, które to robią (wszystko w Pythonie byłoby świetne).

Moja próba

Poniżej znajduje się opis tego, czego próbowałem i problemów, które zauważyłem. Jest to trochę długa lektura i przepraszam, jeśli zawiera warunki, które zostały niewłaściwie użyte. Jeśli można wskazać, co jest błędne, byłoby to bardzo pomocne. Ale moje pytanie jest pogrubione powyżej.

1. PSD i CSD można zapisać jako wartość oczekiwaną (lub średnią zespoloną) produktów przekształceń Fouriera szeregów czasowych. Tak więc macierz kowariancji można zapisać jako: gdzie $$\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. Macierz kowariancji jest macierzą hermitowską, mającą rzeczywiste wartości własne, które są zerowe lub dodatnie. Można go więc rozłożyć na gdzie jest macierzą diagonalną, której niezerowe elementy są pierwiastkami kwadratowymi wartości własnych \ mathbf {C} (f) ; \ mathbf {X} (f) jest macierzą, której kolumny są ortonormalnymi wektorami własnymi \ mathbf {C} (f) ; $$\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}$ to macierz tożsamości.
3. Macierz tożsamości jest zapisywana jako gdzie i to nieskorelowane i złożone szeregi częstotliwości o zerowej średniej i wariancji jednostkowej. $$\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. Używając 3. w 2., a następnie porównaj z 1. Transformaty Fouriera szeregów czasowych to: $$\mathbf{Y}(f) = \sqrt{ \frac{\tau}{2} } \mathbf{z}(f) \: \boldsymbol\lambda^{\frac{1}{2}}(f) \: \mathbf{X}^{\dagger}(f) \;.$$
5. Szeregi czasowe można następnie uzyskać za pomocą procedur takich jak odwrotna szybka transformata Fouriera.

Napisałem procedurę w Pythonie, aby to zrobić:

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

Zastosowałem tę procedurę do PSD i CSD, których analityczne wyrażenia zostały uzyskane z modelowania jakiegoś detektora, z którym pracuję. Ważną rzeczą jest to, że przy wszystkich częstotliwościach tworzą one macierz kowariancji (przynajmniej przynajmniej przekazują wszystkie te ifrutynowe czynności). Macierz kowariancji wynosi 3x3. 3 szeregi czasowe wygenerowano około 9000 razy, a szacowane PSD i CSD, uśrednione dla wszystkich tych realizacji, są przedstawione poniżej z analitycznymi. Chociaż ogólne kształty są zgodne, w CDPW zauważalne są hałaśliwe funkcje przy określonych częstotliwościach (ryc. 2). Po zbliżeniu wokół szczytów na PSD (ryc. 3) zauważyłem, że PSD są w rzeczywistości niedocenianeoraz że hałaśliwe funkcje w dyskach CSD występują prawie w przybliżeniu na tych samych częstotliwościach, co szczyty w dyskach PSD. Nie sądzę, że to zbieg okoliczności i że w jakiś sposób moc wycieka z PSD do CSD. Spodziewałbym się, że krzywe leżą jeden na drugim, przy tak wielu realizacjach danych.

Ponieważ sygnały są nieruchome, prostym podejściem byłoby wykorzystanie białego szumu jako podstawy i przefiltrowanie go, aby pasował do PSD. Sposobem na obliczenie tych współczynników filtrowania jest zastosowanie prognozy liniowej .

Wygląda na to, że istnieje dla niego funkcja python, wypróbuj:

from scikits.talkbox import lpc

Jeśli chcesz (użyłem tylko odpowiednika MATLAB). Jest to podejście stosowane w przetwarzaniu mowy, w którym formanty są szacowane w ten sposób.

Trochę za późno na imprezę, jak zwykle, ale widzę trochę niedawnej aktywności, więc będę miał dwa jeny.

Po pierwsze, nie mogę winić próby OP - wygląda mi to dobrze. Rozbieżności mogą być spowodowane problemami ze skończonymi próbkami, na przykład dodatnim błędem oszacowania mocy sygnału.

Myślę jednak, że istnieją prostsze sposoby generowania szeregów czasowych z macierzy gęstości widma krzyżowego (CPSD, to właśnie OP nazywa macierzą kowariancji).

Jednym podejściem parametrycznym jest użycie CPSD w celu uzyskania opisu autoregresji, a następnie użycie go do wygenerowania szeregów czasowych. W Matlab możesz to zrobić za pomocą narzędzi przyczynowości Grangera (np. Multivaraite Granger Causality Toolbox, Seth, Barnett ). Zestaw narzędzi jest bardzo łatwy w użyciu. Ponieważ istnienie CPSD gwarantuje autoregresyjny opis, podejście to jest dokładne. (aby uzyskać więcej informacji na temat CPSD i autoregresji, patrz „Pomiar zależności liniowej i sprzężenie zwrotne między szeregami czasowymi” Geweke, 1982 lub wiele artykułów Aneil Seth + Lionel Barnett, aby uzyskać pełny obraz).

Potencjalnie prościej jest zauważyć, że CPSD można utworzyć, stosując fft do automatycznej kowariancji (podając przekątną CPSD, tj. Moc sygnałów) i kowariancję krzyżową (podając elementy przekątne, tj. Moc krzyżową). Zatem poprzez zastosowanie odwrotnego fft do CPSD możemy uzyskać autokorelację i automatyczną kowariancję. Następnie możemy je wykorzystać do wygenerowania próbek naszych danych.

Mam nadzieję że to pomoże. Proszę zostawić wszelkie prośby o informacje w komentarzach, a ja postaram się odpowiedzieć.