Как взять преобразование Фурье от дискретных точек?

Question

[Сергей Соколов]

Временные данные представляют из себя просто отсчёты времени, моменты, в которые прозошло событие. Напр. события отмечены в 1-ю секунду, 3, 5 и т.д.: t = [1, 3, 5, 7, 9].

Как правильно получить частотный спектр этих данных? Ожидаю, что должен вылезти пик на амплитуде "2".

Пробовал выбрать некую частоту дискретизации, напр. с шагом в 0.5 секунды, и сделать расчёску из "0", а там, где были события – "1". И от этого набора брать fft()

Но получается совсем не то, что ожидал:


Мой код Matlab этой попытки:

A=[1,3,5,7,9];
subplot(2,1,1);
B = zeros(max(A),1);
B(A) = 1;
bar(B);

dx = 0.5;
X = zeros( round( max(A) / dx) - 1, 1);
X( round(A / dx)) = 1;
t = 1:length(X);
Y = fft(X);
P2 = abs(Y/length(X));
L = length(X);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = 0:L/2;
subplot(2,1,2);
plot(f,P1,'o-');
axis([0 10 0 0.5]);

Answer 1

[longclaps]

> Как правильно получить частотный спектр этих данных?
Спектр всегда частотный, а масло - масляное
> Ожидаю, что должен вылезти пик на амплитуде "2".
Вот это вообще непонятно что означает

Смотри: у тебя есть какие-то дискретные события, происходящие в произвольные моменты времени, например, ты сидишь под яблоней и, когда падает яблоко - засекаешь время. И ты хочешь что-то понять про периодичность, связанную с этими событиями.
Первое - временные рамки. Помимо времени падения первого и последнего яблока понадобятся моменты начала и конца наблюдений. Пусть это будут 00:00 и 24:00.
Второе. Дискретное преобразование фурье работает с периодическими отсчетами, например, сколько яблок упало за каждый час в течение суток (те 24 отсчета). БПФ - алгоритм, требующий, чтобы число отсчетов было степенью двойки (в моём примере возьму 32 отсчета по 45 мин).
Третье. Если яблок падает мало, а вытянуть из опыта хочется по-максимуму, можно попробовать интерполяцию, например, если в 12:30 упало 1 яблоко, это можно интерпретировать так: в 12:00 упала 1/3 яблока и в 12:45 - остальные 2/3. Для корректной интерполяции на краях временного диапазона неплохо бы прихватить по 45мин наблюдений до и после.

Answer 2

[nirvimel]

Не подскажу про Matlab, но на Python это делается так:

import numpy
import scipy.fftpack
x = numpy.array([0, 1, 0, 1, 0, 1, 0, 1, 0, 1])
numpy.abs(scipy.fftpack.fft(x)[:(len(x)/2)]) * (2./len(x))

На выходе: array([ 1., 0., 0., 0., 0.]) - тут виден единственный пик на самой высокой частоте n все более низкие частоты по нулям.
Если взять другой пример:

x = numpy.array([0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1])

, то получим array([ 0.5, 0. , 0. , 0.5, 0. , 0. ]) - тут два равных пика на n и на n/4, последний очевидно следует из данных, первый из того факта, что фигура исходного сигнала далека от синусоиды.

Answer 3

[Сергей]

Я чего-то давно не практиковал Matlab. Лет 10 :) Все забыл.
Но скажу так -- преобразование фурье должно давать вам комплексное число, соотвественно пик вам нужно наблюдать на АЧХ (модуль комплексного числа).

Потом ваш сигнал - это меандр.
Может это поможет: https://habrahabr.ru/post/219337/