Mensimulasikan deret waktu yang diberikan daya dan lintas kerapatan spektral

Saya mengalami masalah dalam menghasilkan serangkaian seri waktu berwarna stasioner, mengingat matriks kovariansnya (densitas spektral daya (PSD)) dan densitas spektral cross-power (CSD)).

Saya tahu bahwa, mengingat dua seri waktu $y_{I}(t)$ dan $y_{J}(t)$ , saya dapat memperkirakan kepadatan spektral daya (PSD) dan kepadatan spektral silang (CSD) menggunakan banyak rutinitas yang tersedia secara luas, seperti psd()dan csd()fungsi dalam Matlab , dll. PSD dan CSD membentuk matriks kovarian:

$$\mathbf{C}(f) = \left( \begin{array}{cc} P_{II}(f) & P_{IJ}(f)\\ P_{JI}(f) & P_{JJ}(f) \end{array} \right)\;,$$ yang pada umumnya merupakan fungsi frekuensi$f$ .

Apa yang terjadi jika saya ingin melakukan yang sebaliknya? Mengingat matriks kovarians, bagaimana cara menghasilkan realisasi $y_{I}(t)$ dan $y_{J}(t)$ ?

Harap sertakan teori latar belakang apa pun, atau tunjukkan alat yang ada yang melakukan ini (apa pun dengan Python akan bagus).

Usaha saya

Di bawah ini adalah deskripsi dari apa yang saya coba, dan masalah yang saya perhatikan. Ini agak lama dibaca, dan maaf jika mengandung istilah yang telah disalahgunakan. Jika apa yang salah dapat ditunjukkan, itu akan sangat membantu. Tetapi pertanyaan saya adalah yang dicetak tebal di atas.

1. PSD dan CSD dapat ditulis sebagai nilai ekspektasi (atau rata-rata ensemble) dari produk-produk dari transformasi Fourier dari deret waktu. Jadi, matriks kovarians dapat ditulis sebagai: $$\mathbf{C}(f) = \frac{2}{\tau} \langle \mathbf{Y}^{\dagger}(f) \mathbf{Y}(f) \rangle \;,$$ mana $$\mathbf{Y}(f) = \left( \begin{array}{cc} \tilde{y}_{I}(f) & \tilde{y}_{J}(f) \end{array} \right) \;.$$
2. Matriks kovarians adalah matriks Hermitian, yang memiliki nilai eigen nyata yang nol atau positif. Jadi, ini dapat didekomposisi menjadi mana adalah matriks diagonal yang unsur-unsurnya bukan nol adalah akar kuadrat dari nilai eigen ; adalah matriks yang kolomnya adalah vektor eigen ortonormal dari ; $$\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}$ adalah matriks identitas.
3. Matriks identitas ditulis sebagai mana dan tidak berkorelasi dan seri-frekuensi kompleks dengan nol rata-rata dan varians unit. $$\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. Dengan menggunakan 3. pada 2., dan kemudian membandingkan dengan 1. Transformasi Fourier dari deret waktu adalah: $$\mathbf{Y}(f) = \sqrt{ \frac{\tau}{2} } \mathbf{z}(f) \: \boldsymbol\lambda^{\frac{1}{2}}(f) \: \mathbf{X}^{\dagger}(f) \;.$$
5. Rangkaian waktu kemudian dapat diperoleh dengan menggunakan rutinitas seperti transformasi Fourier cepat terbalik.

Saya telah menulis rutin dalam Python untuk melakukan ini:

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

Saya telah menerapkan rutin ini untuk PSD dan CSD, ekspresi analitik yang telah diperoleh dari pemodelan beberapa detektor yang saya kerjakan. Yang penting adalah bahwa pada semua frekuensi, mereka membentuk matriks kovarians (setidaknya mereka melewati semua ifpernyataan itu dalam rutinitas). Matriks kovarians adalah 3x3. Seri-3 waktu telah dihasilkan sekitar 9000 kali, dan perkiraan PSD dan CSD, dirata-rata untuk semua realisasi ini diplot di bawah ini dengan yang analitis. Sementara bentuk keseluruhan setuju, ada fitur berisik yang terlihat pada frekuensi tertentu di CSD (Gbr.2). Setelah close-up di sekitar puncak di PSD (Gbr.3), saya perhatikan bahwa PSD sebenarnya diremehkan, dan bahwa fitur bising di CSD terjadi pada frekuensi yang hampir sama dengan puncak di PSD. Saya tidak berpikir bahwa ini adalah kebetulan, dan entah bagaimana kekuatan bocor dari PSD ke CSD. Saya berharap kurva berada di atas satu sama lain, dengan banyak realisasi data.

Karena sinyal Anda diam, pendekatan sederhana adalah menggunakan white noise sebagai dasar dan memfilternya agar sesuai dengan PSD Anda. Cara untuk menghitung koefisien filter ini adalah dengan menggunakan prediksi linier .

Tampaknya ada fungsi python untuk itu, coba saja:

from scikits.talkbox import lpc

Jika Anda mau (saya hanya menggunakan setara MATLAB). Ini adalah pendekatan yang digunakan dalam pemrosesan wicara, di mana forman diperkirakan dengan cara ini.

Agak terlambat ke pesta, seperti biasa, tapi aku melihat beberapa kegiatan baru-baru ini jadi aku akan dua yen ku

Pertama, saya tidak bisa menyalahkan upaya OP - itu terlihat benar bagi saya. Perbedaan ini dapat disebabkan oleh masalah dengan sampel yang terbatas, misalnya bias positif estimasi daya sinyal.

Namun, saya berpikir bahwa ada cara yang lebih sederhana untuk menghasilkan deret waktu dari matriks kerapatan spektral silang (CPSD, inilah yang disebut OP sebagai matriks kovarians).

Salah satu pendekatan parametrik adalah menggunakan CPSD untuk mendapatkan deskripsi autoregresif dan kemudian menggunakannya untuk menghasilkan deret waktu. Di matlab Anda dapat melakukan ini dengan menggunakan alat kausalitas Granger (mis. Toolalitas kausalitas Granger Multivaraite, Seth, Barnett ). Toolbox sangat mudah digunakan. Karena keberadaan CPSD menjamin deskripsi autoregresif pendekatan ini tepat. (untuk info lebih lanjut tentang CPSD dan autoregresi lihat "Pengukuran Ketergantungan Linier dan Umpan Balik antara Berbagai Waktu" oleh Geweke, 1982, atau banyak makalah Aneil Seth + Lionel Barnett, untuk mendapatkan gambaran lengkap).

Berpotensi lebih sederhana, mencatat CPSD dapat dibentuk dengan menerapkan fft ke kovarians otomatis (memberikan diagonal CPSD, yaitu kekuatan sinyal) dan kovarians silang (memberikan elemen diagonal off, yaitu kekuatan lintas). Jadi dengan menerapkan invers fft ke CPSD kita bisa mendapatkan autokorelasi dan kovarian otomatis. Kami kemudian dapat menggunakan ini untuk menghasilkan sampel data kami.

Semoga ini membantu. Silakan tinggalkan permintaan info di komentar dan saya akan mencoba untuk mengatasinya.