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:
psd()csd()
Co się stanie, jeśli chcę zrobić odwrotnie? Biorąc pod uwagę macierz kowariancji, jak wygenerować realizację i ?
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.
- 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
- 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) ;to macierz tożsamości.
- 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.
- Używając 3. w 2., a następnie porównaj z 1. Transformaty Fouriera szeregów czasowych to:
- 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.
