Необъяснимая симметрия при вычислении спектральной плотности мощности белого шума

Я пытаюсь узнать больше о шуме, спектральной плотности мощности (PSD) и статистических отклонениях. Что касается этого, я пытаюсь вычислить спектральную плотность мощности белого шума, однако, когда я это делаю, я получаю очень странную симметрию. мой спектр кажется симметричным относительно значения центральной частоты, что явно неверно. Я новичок в использовании автокорреляций и спектральной плотности мощности, поэтому я был бы признателен, если бы кто-нибудь помог мне подтолкнуть меня к ошибке.

Код для расчета PSD:

import numpy as np 
from math import sin, pi, log10
from allan_noise import white,pink,brown, violet
import acor
import numpy as np


#METHOD ONE: AUTOCORRELATION FUNCTION METHOD
def psd1(time_dats):
    #auto_corr = acor.function(time_dats)
    auto_corr = np.correlate(time_dats,time_dats,mode="same")

    #To check autocorrelation
    for i in range(len(auto_corr)):
        print i*.0000001 ,auto_corr[i]

    return fft(auto_corr) 

#DEFINE VARIABLES
t = .0001       #simulation time
N = 100000  #number of data points 
dt = t/N    #time interval between data points

#CREATE SIGNAL
sig = white(N)
df = 1.0/(t)
freq_axis = np.arange(0,N/t,df)

spec = psd1(sig)

#OPEN UP OUTPUT FILE
f = open('data/psdtest_f','w')
g = open('data/psdtest_t','w')

#PRINT OUT DATA
for i in range(N):
    f.write('%g %g\n' %(freq_axis[i],log10(abs(spec[i]))))
    g.write('%g %g\n' %(i*dt, sig[i]))

Используя этот код, я создаю следующие графики, доступ к которым можно получить здесь https://drive.google.com/#folders/0B8BdiIhqTYw7Vk1EaDFMQW84RHM :

  1. Временной профиль шума до расчетов

  2. Автокорреляционная функция, рассчитанная по временному профилю (я знаю, что шкала по оси X неверна, но это не влияет на код где-либо еще

  3. Power Spectral Denisty, красиво симметричный, хотя и не должен быть

Любая помощь, которую кто-либо мог бы оказать относительно того, что вызывает эту симметрию, была бы наиболее полезной. я протестировал код с простым синусоидальным сигналом и получил ожидаемые результаты (без симметрии)


python signals symmetric
person Leavenotrace    schedule 19.06.2014    source источник
comment
Этот вопрос кажется не по теме, потому что он касается интерпретации результатов, а не реализации кода.   -  person jonrsharpe    schedule 19.06.2014
comment
Это может быть лучше подходит, например. Обработка сигналов или, возможно, Перекрестная проверка   -  person jonrsharpe    schedule 19.06.2014
comment
Извиняюсь, хотя я чувствую, что в реализации кода я ошибаюсь. Я не совсем понимаю, как работают функции numpy. хотя, если вы чувствуете, что это не по теме, я, конечно, могу перенести это   -  person Leavenotrace    schedule 19.06.2014
comment
Вы считаете, что реализация неверна? Каких результатов вы ожидали? Вы тестировали код - что работает, что нет? Запрос на объяснение кода будет здесь не по теме, а на общий вопрос что происходит не так? сложно ответить — см. Как спросить и рассмотрите возможность предоставления минимального примера   -  person jonrsharpe    schedule 19.06.2014
comment
Вы правы, я новичок в stackoverflow и все еще привыкаю к ​​нему. Я не прошу никого объяснять, как работает библиотека numpy, я просто говорю, что думаю, что это может быть моя ошибка, и кто-то с немного большим опытом может ее подобрать. Я попытался указать внизу, что я протестировал код с использованием функции sin и получил пик, соответствующий одной частоте. Белый шум должен давать плоскую спектральную плотность мощности. в основном это должно выглядеть как первая половина графика, который я предоставил, что наводит меня на мысль, что я где-то отражаю данные.   -  person Leavenotrace    schedule 19.06.2014
comment
Вы не получаете какой-то граничный эффект? Если вы сигнализируете, и то, что вы коррелируете, не полностью перекрывается, у вас будут меньшие значения на краях, поскольку в них вносится меньше значений.   -  person Salix alba    schedule 19.06.2014

Ответы (1)


arrow_upward
0
arrow_downward

Во-первых, ваша функция написана неправильно, вы берете fft из автозамены, которая не является вашим исходным массивом сигналов. Забавно, что из-за характера ошибки округления вы также получите шум, такой как psd. Во-вторых, вы неправильно вычисляете ось частот, так как они должны быть расширены до N/(t*2) (это фр. Найквиста). Вместо этого, поскольку длина вашего массива freq_axis равна N, вы извлекаете N элементов из своего массива sig, и из-за этого вы просто читаете отрицательные частоты с тем же psd, что вызывает у вас симметрию (взяв журнал, преобразующий вашу переменную в реальную и все Фурье коэффициенты для отрицательных частот просто сопряжены для положительных). замена вашего кода на следующий дает отличный результат:

sig = white(N,1,N/t)
(siglog,freq_axis)=ml.psd(sig,N,(N/t), detrend=ml.detrend_none,
   window=ml.window_hanning, noverlap=0, pad_to=None,
    sides='default', scale_by_freq=None)
plt.plot(freq_axis,np.log10(siglog))
plt.show()

результат matplotlib.mlab.psd

не забудьте импортировать

import matplotlib.mlab as ml
person Andrey Debolskiy    schedule 29.01.2016