Какой самый быстрый способ приблизить период данных, используя Octave? - PullRequest
Новая
       24

Какой самый быстрый способ приблизить период данных, используя Octave?

4 голосов
/ 03 апреля 2010

У меня есть набор данных, который является периодическим (но не синусоидальным). У меня есть набор значений времени в одном векторе и набор амплитуд во втором векторе. Я хотел бы быстро приблизить период функции. Есть предложения?

В частности, вот мой текущий код. Я хотел бы приблизить период вектора x (:, 2) к вектору t. В конечном счете, я хотел бы сделать это для множества начальных условий, рассчитать период каждого и построить график результата. 

function xdot = f (x,t)
         xdot(1) =x(2);
         xdot(2) =-sin(x(1));
endfunction

x0=[1;1.75];     #eventually, I'd like to try lots of values for x0(2)
t = linspace (0, 50, 200);


x = lsode ("f", x0, t)

plot(x(:,1),x(:,2));

Спасибо!

John

Ответы [ 2 ]

6 голосов
/ 04 апреля 2010

Посмотрите на функцию автокорреляции.

Из Википедия

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

Пол Бурк описывает, как эффективно рассчитать функцию автокорреляции на основе быстрого преобразования Фурье ( link ).

2 голосов
/ 04 апреля 2010

Дискретное преобразование Фурье может дать вам периодичность. Более длинное временное окно дает вам большее разрешение по частоте, поэтому я изменил ваше определение t на t = linspace(0, 500, 2000). временной домен http://img402.imageshack.us/img402/8775/timedomain.png (вот ссылка на сюжет , на хостинг-сайте она выглядит лучше). Вы могли бы сделать: 

h = hann(length(x), 'periodic'); %# use a Hann window to reduce leakage
y = fft(x .* [h h]); %# window each time signal and calculate FFT
df = 1/t(end); %# if t is in seconds, df is in Hz
ym = abs(y(1:(length(y)/2), :)); %# we just want amplitude of 0..pi frequency components
semilogy(((1:length(ym))-1)*df, ym);

частотная область http://img406.imageshack.us/img406/2696/freqdomain.png Участок участка.

Глядя на график, первый пик находится на отметке 0,06 Гц, что соответствует 16-секундному периоду, замеченному в plot(t,x).

Хотя это не так быстро в вычислительном отношении. БПФ - это N * log (N) операций.

...