Непрерывное преобразование Фурье на дискретных данных с использованием Mathematica?

Question

У меня есть некоторые периодические данные, но объем данных не кратен периоду.Как я могу Фурье проанализировать эти данные?Пример:

% Давайте создадим некоторые данные для тестирования:

data = Table[N[753+919*Sin[x/623-125]], {x,1,25000}] 

% Я сейчас получаю эти данные, но не знаю, что они получены из формулы выше.Я пытаюсь восстановить формулу только из «данных».

% Рассматривая первые несколько непостоянных членов ряда Фурье:

ListPlot[Table[Abs[Fourier[data]][[x]], {x,2,20}], PlotJoined->True, 
 PlotRange->All] 

показывает ожидаемый всплеск в 6 (поскольку числопериоды действительно 25000 / (623 * 2 * Pi) или около 6,38663, хотя мы этого не знаем).

% Теперь, как мне вернуться к 6.38663?Одним из способов является «свертка» данных с произвольными кратными Cos [x].

convolve[n_] := Sum[data[[x]]*Cos[n*x], {x,1,25000}] 

% И изобразите «свертку» вблизи n = 6:

Plot[convolve[n],{n,5,7}, PlotRange->All] 

, мы видим всплеск примерно там, где и ожидалось.

% Мы пытаемся найти FindMaximum:

FindMaximum[convolve[n],{n,5,7}] 

, но результат бесполезен и неточен:

FindMaximum::fmmp:  
   Machine precision is insufficient to achieve the requested accuracy or 
    precision. 

Out[119]= {98.9285, {n -> 5.17881}} 

, потому что функция очень волнистая.

% Уточняя наш интервал (используя визуальный анализ на графиках), мы, наконец, находим интервал, в котором свертка [] не слишком сильно покачивается:

Plot[convolve[n],{n,6.2831,6.2833}, PlotRange->All] 

и FindMaximum работают:

FindMaximum[convolve[n],{n,6.2831,6.2833}] // FortranForm 
List(1.984759605826571e7,List(Rule(n,6.2831853071787975))) 

% Тем не менее, этот процесс уродлив, требует вмешательства человека, и вычисление convolve [] ДЕЙСТВИТЕЛЬНО медленное.Есть лучший способ сделать это?

% Глядя на ряд данных Фурье, могу ли я как-то угадать "истинное" количество периодов, равное 6,38663?Конечно, фактический результат будет 6,283185, так как мои данные соответствуют этому лучше (потому что я выбираю только в конечном количестве точек).

Answer 1

На основе справки Mathematica для функции Фурье / Приложения / Идентификация частоты: Проверено на версии 7

n = 25000;
data = Table[N[753 + 919*Sin[x/623 - 125]], {x, 1, n}];
pdata = data - Total[data]/Length[data];
f = Abs[Fourier[pdata]];
pos = Ordering[-f, 1][[1]]; (*the position of the first Maximal value*)  
fr = Abs[Fourier[pdata Exp[2 Pi I (pos - 2) N[Range[0, n - 1]]/n], 
   FourierParameters -> {0, 2/n}]];
frpos = Ordering[-fr, 1][[1]];

N[(pos - 2 + 2 (frpos - 1)/n)]

возвращает 6,37072

Answer 2

Ищите длину периода, используя автокорреляцию, чтобы получить оценку:

autocorrelate[data_, d_] := 
 Plus @@ (Drop[data, d]*Drop[data, -d])/(Length[data] - d)

ListPlot[Table[{d, autocorrelate[data, d]}, {d, 0, 5000, 100}]]

Интеллектуальный поиск первого максимума вдали от d = 0 может быть наилучшей оценкой, которую вы можете получить из доступных данных?

Answer 3

(* the data *) 

data = Table[N[753+919*Sin[x/623-125]], {x,1,25000}]; 

(* Find the position of the largest Fourier coefficient, after 
removing the last half of the list (which is redundant) and the 
constant term; the [[1]] is necessary because Ordering returns a list *) 

f2 = Ordering[Abs[Take[Fourier[data], {2,Round[Length[data]/2+1]}]],-1][[1]] 

(* Result: 6 *) 

(* Directly find the least squares difference between all functions of 
the form a+b*Sin[c*n-d], with intelligent starting values *) 

sol = FindMinimum[Sum[((a+b*Sin[c*n-d]) - data[[n]])^2, {n,1,Length[data]}], 
{{a,Mean[data]},{b,(Max[data]-Min[data])/2},{c,2*f2*Pi/Length[data]},d}] 

(* Result (using //InputForm):  

FindMinimum::sszero:  
   The step size in the search has become less than the tolerance prescribed by 
   the PrecisionGoal option, but the gradient is larger than the tolerance 
   specified by the AccuracyGoal option. There is a possibility that the method 
   has stalled at a point that is not a local minimum. 

{2.1375902350021628*^-19, {a -> 753., b -> -919., c -> 0.0016051364365971107,  
  d -> 2.477886509998064}} 

*) 


(* Create a table of values for the resulting function to compare to 'data' *) 

tab = Table[a+b*Sin[c*x-d], {x,1,Length[data]}] /. sol[[2]]; 

(* The maximal difference is effectively 0 *) 

Max[Abs[data-tab]] // InputForm 

(* Result: 7.73070496506989*^-12 *) 

Хотя вышесказанное не обязательно полностью отвечает на мой вопрос, я нашел его несколько примечательным.

Раньше я пытался использовать FindFit[] с Method -> NMinimize (что должно обеспечить лучшее глобальное соответствие), но это не сработало, возможно потому, что вы не можете дать FindFit[]интеллектуальные стартовые значения.

Ошибка, которую я получаю, беспокоит меня, но, похоже, не имеет значения.