Моделирование временных рядов с учетом мощности и кросс-спектральных плотностей

20

У меня возникают проблемы при создании набора стационарных цветных временных рядов, учитывая их ковариационную матрицу (их спектральные плотности мощности (PSD) и спектральные плотности перекрестных мощностей (CSD)).

Я знаю, что, учитывая два временных ряда и , я могу оценить их спектральные плотности мощности (PSD) и кросс-спектральные плотности (CSD), используя многие широко доступные процедуры, такие как и функции в Matlab и т. д. PSD и CSD составляют ковариационную матрицу: yI(t)yJ(t)psd()csd()

С(е)знак равно(пяя(е)пяJ(е)пJя(е)пJJ(е)),
который в общем случае является функцией частоты . е

Что произойдет, если я захочу сделать наоборот? Учитывая ковариационную матрицу, как мне сгенерировать реализацию и ?Yя(T)YJ(T)

Пожалуйста, включите любую фоновую теорию или укажите какие-либо существующие инструменты, которые делают это (все в Python было бы здорово).

Моя попытка

Ниже приведено описание того, что я пробовал, и проблем, которые я заметил. Это немного долго читать, и извините, если он содержит термины, которые были использованы неправильно. Если можно указать на то, что является ошибочным, это было бы очень полезно. Но мой вопрос выделен жирным шрифтом выше.

  1. PSD и CSD могут быть записаны как ожидаемое значение (или среднее по ансамблю) произведений преобразований Фурье временного ряда. Таким образом, ковариационная матрица может быть записана как: где
    С(е)знак равно2τY(е)Y(е),
    Y(е)знак равно(Y~я(е)Y~J(е)),
  2. Ковариационная матрица - это эрмитова матрица, имеющая действительные собственные значения, которые либо равны нулю, либо положительны. Таким образом, его можно разложить на где \ lambda ^ {\ frac {1} {2}} ( f) - диагональная матрица, ненулевыми элементами которой являются квадратные корни из собственных значений \ mathbf {C} (f) ; \ mathbf {X} (f) - это матрица, столбцы которой являются ортонормированными собственными векторами \ mathbf {C} (f) ;
    С(е)знак равноИкс(е)λ12(е)яλ12(е)Икс(е),
    λ12(е)С(е)Икс(е)С(е)я - это единичная матрица.
  3. Тождественная матрица записывается как где и - некоррелированные и сложные частотные ряды с нулевым средним и единичной дисперсией.
    язнак равноZ(е)Z(е),
    Z(е)знак равно(Zя(е)ZJ(е)),
    {Zя(е)}язнак равноя,J
  4. Используя 3. в 2., а затем сравните с 1. Преобразования Фурье временного ряда:
    Y(е)знак равноτ2Z(е)λ12(е)Икс(е),
  5. Временные ряды могут затем быть получены с использованием процедур, таких как обратное быстрое преобразование Фурье.

Я написал процедуру на Python для этого:

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

Я применил эту процедуру к PSD и CSD, аналитические выражения которых были получены при моделировании детектора, с которым я работаю. Важно то, что на всех частотах они составляют ковариационную матрицу (по крайней мере, они пропускают все эти ifутверждения в процедуре). Ковариационная матрица 3х3. 3 временных ряда были сгенерированы около 9000 раз, и оценочные PSD и CSD, усредненные по всем этим реализациям, представлены ниже с аналитическими. В то время как общие формы согласуются, на определенных частотах в CSD есть заметные шумные особенности (Рис.2). После изучения пиков в PSD (рис. 3) я заметил, что PSD фактически недооцененыи что шумовые характеристики в CSD возникают примерно на тех же частотах, что и пики в PSD. Я не думаю, что это совпадение, и что каким-то образом энергия просачивается из PSD в CSD. Я ожидал бы, что кривые будут лежать друг на друге с таким большим количеством реализаций данных.

Рисунок 1: P11
Рисунок 2: P12 Рисунок 2: P11 (крупный план)

я женат
источник
Добро пожаловать на сайт. Я частично проголосовал за этот вопрос, чтобы вы не могли публиковать изображения. Если нет, просто опубликуйте ссылки, и кто-то с достаточной репутацией отредактирует для вставки изображений.
кардинал
1
Вы пытались отфильтровать высокочастотный шум?
Карл

Ответы:

1

Поскольку ваши сигналы являются стационарными, простой подход состоит в том, чтобы использовать белый шум в качестве основы и фильтровать его для соответствия вашим PSD. Способ вычисления этих коэффициентов фильтра состоит в использовании линейного прогнозирования .

Кажется, для этого есть функция python, попробуйте:

from scikits.talkbox import lpc

Если хотите (я использовал только эквивалент MATLAB). Этот подход используется в обработке речи, где форманты оцениваются таким образом.

Йонас Шварц
источник
Разве вы не имеете в виду применять фильтр к сигналу, а не к белому шуму?
Майкл Р. Черник
Нет, я стремлюсь приблизиться к фильтру, в котором передаточная функция напоминает PSD стационарного процесса. Если белый шум, который имеет одинаковую мощность во всех полосах частот, фильтруется с ними, выходной сигнал оптимально напоминает исходные сигналы по спектральной плотности мощности.
Йонас Шварц
0

Немного опоздал на вечеринку, как обычно, но я вижу некоторую недавнюю активность, поэтому я получу свои две иены.

Во-первых, я не могу винить в попытках ОП - мне это кажется правильным. Расхождения могут быть вызваны проблемами с конечными выборками, например, положительным смещением оценки мощности сигнала.

Тем не менее, я думаю, что есть более простые способы генерирования временных рядов из матрицы кросс-спектральной плотности (CPSD, это то, что OP называет ковариационной матрицей).

Одним из параметрических подходов является использование CPSD для получения авторегрессионного описания, а затем его использование для генерации временных рядов. В Matlab вы можете сделать это, используя инструменты причинности Грейнджер (например, набор инструментов причинности Мультиварайт Грейнджер, Сет, Барнетт ). Панель инструментов очень проста в использовании. Поскольку существование CPSD гарантирует авторегрессионное описание, этот подход является точным. (для получения дополнительной информации о CPSD и авторегрессии см. «Измерение линейной зависимости и обратной связи между несколькими временными рядами», Geweke, 1982, или многие из работ Aneil Seth + Lionel Barnett, чтобы получить полную картину).

Потенциально проще отметить, что CPSD может быть сформирован путем применения fft к автоковариантности (давая диагональ CPSD, то есть мощность сигналов) и кросс-ковариации (давая недиагональные элементы, то есть перекрестную мощность). Таким образом, применяя обратное fft к CPSD, мы можем получить автокорреляцию и автоковариацию. Затем мы можем использовать их для создания образцов наших данных.

Надеюсь это поможет. Пожалуйста, оставляйте любые запросы на информацию в комментариях, и я постараюсь обратиться.

dcneuro
источник