Есть ли алгоритм для вычисления фазы для одной частоты?

17

Если у вас есть функция и эталонная волна что будет быстрым алгоритмом для вычисления ?f(t)=Asin(ωt+ϕ)sin(ωx)ϕ

Я смотрел на алгоритм Гёртцела , но он не имеет отношения к фазе?

SamFisher83
источник

Ответы:

5

Используйте ДПФ на определенной частоте. Затем вычислите амплитуду и фазу из реальной / воображаемой частей. Это дает вам фазу, относящуюся к началу времени выборки.

В «нормальном» БПФ (или ДПФ, рассчитанном для всех N гармоник) вы обычно вычисляете частоту с помощью f = k * (sample_rate) / N, где k - целое число. Хотя это может показаться кощунственным (особенно членам Церкви Полностью Целого), вы можете использовать нецелые значения k при выполнении одного ДПФ.

Например, предположим, что вы сгенерировали (или получили) N = 256 точек синусоиды 27 Гц. (скажем, sample_rate = 200). Ваши «нормальные» частоты для 256-точечного БПФ (или N-точечного БПФ) будут соответствовать: f = k * (sample_rate) / N = k * (200) / 256, где k - целое число. Но нецелое «k» 34,56 будет соответствовать частоте 27 Гц, используя параметры, перечисленные выше. Это похоже на создание «корзины» DFT, которая точно центрирована на интересующей частоте (27 Гц.). Некоторый код C ++ (DevC ++ компилятор) может выглядеть следующим образом:

#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <cmath>
using namespace std;

// arguments in main needed for Dev-C++ I/O
int main (int nNumberofArgs, char* pszArgs[ ] ) {
const long N = 256 ;
double sample_rate = 200., amp, phase, t, C, S, twopi = 6.2831853071795865; 
double  r[N] = {0.}, i[N] = {0.}, R = 0., I = 0. ;
long n ;

// k need not be integer
double k = 34.56;

// generate real points
for (n = 0; n < N; n++) {
    t =  n/sample_rate;
    r[n] = 10.*cos(twopi*27.*t - twopi/4.);
}  // end for

// compute one DFT
for (n = 0; n < N; n++) {
    C = cos(twopi*n*k/N); S = sin(twopi*n*k/N);
    R = R + r[n]*C + i[n]*S;
    I = I + i[n]*C - r[n]*S;
} // end for

cout<<"\n\ndft results for N = " << N << "\n";
cout<<"\nindex k     real          imaginary       amplitude         phase\n";

amp = 2*sqrt( (R/N)*(R/N) + (I/N)*(I/N) ) ;
phase = atan2( I, R ) ;
// printed R and I are scaled
printf("%4.2f\t%11.8f\t%11.8f\t%11.8f\t%11.8f\n",k,R/N,I/N,amp,phase);

cout << "\n\n";
system ("PAUSE");
return 0;
} // end main

//**** end program

(PS: я надеюсь, что вышеупомянутое хорошо транслируется в stackoverflow - некоторые из них могут обернуться)

Результатом вышеупомянутого является фаза -twopi / 4, как показано в сгенерированных реальных точках (и усиление удваивается для отражения частоты положительного / отрицательного частот).

Несколько вещей, на которые стоит обратить внимание - я использую косинус для генерации тестовой формы волны и интерпретации результатов - вы должны быть осторожны с этим - фаза ссылается на время = 0, то есть когда вы начали сэмплирование (т.е. когда вы собрали r [0] ), а косинус - правильная интерпретация).

Приведенный выше код не является ни элегантным, ни эффективным (например: используйте справочные таблицы для значений sin / cos и т. Д.).

Ваши результаты станут более точными, если вы будете использовать большее N, и есть небольшая ошибка из-за того, что частота дискретизации и N выше не кратны друг другу.

Конечно, если вы хотите изменить частоту дискретизации, N или f, вам придется изменить код и значение k. Вы можете вставить ячейку DFT в любое место на непрерывной частотной линии - просто убедитесь, что вы используете значение k, которое соответствует интересующей частоте.

Кевин МакГи
источник
Этот подход может быть улучшен путем корректировки N, чтобы сделать k ближе к целому. Я разместил отдельный ответ, который ухудшает точность этого алгоритма.
Моджуба
10

Задача может быть сформулирована как (нелинейная) задача наименьших квадратов:

F(ϕ)=12i=1n[Asin(ωi+ϕ)fi(ω)]2

где - целевая функция для минимизации по отношению к .F(ϕ)ϕ

Производная очень проста:

F(ϕ)=i=1nAcos(ωi+ϕ)[Asin(ωi+ϕ)fi(ω)]

Выше целевая функция может быть сведена к минимуму итеративно с помощью метода градиентного спуска (приближение первого порядка), метод Ньютона , метод Гаусса-Ньютон или метод Левенберга-Marquardt (второй порядок аппроксимации - должна быть предоставлены в них).F(ϕ)

Очевидно, что указанная выше целевая функция имеет несколько минимумов из-за периодичности, поэтому можно добавить некоторый штрафной член для различения других минимумов (например, добавив к модельному уравнению). Но я думаю , что оптимизация будет просто сходиться до ближайшего минимума , и вы можете обновить результат вычитания .ϕ22πk,kN

ЛИБОР
источник
Я не думаю, что вам нужно наказывать из-за периодичности нет? Вы можете просто взять любые минимумы в фазовом пространстве, к которому он подходит, и выполнить модуль , нет? 2π
Спейси
@ Mohammad Да, но некоторые методы оптимизации могут использовать несколько начальных точек, которые должны сходиться к одному и тому же значению или предполагать выпуклую функцию с одним глобальным минимизатором, который можно хорошо аппроксимировать квадратичной. Другое преимущество состоит в том, что мы заканчиваем с тем же результатом для любой начальной точки . ϕ0
Либор
Интересный. Могу ли я пригласить вас также ответить на этот вопрос ? :-)
Спейси
@ Мохаммад: Хорошо, я немного внес свой вклад :)
Libor
Куда идет функция fi (w)? fi (w) не является константой, поэтому, когда вы берете производную от непостоянной, как она становится нулевой?
SamFisher83
5

Существует несколько различных формулировок алгоритма Гёртцеля. Те, которые предоставляют 2 переменные состояния (ортогональные или близкие к) или комплексную переменную состояния, в качестве возможных выходных данных часто могут использоваться для вычисления или оценки фазы со ссылкой на некоторую точку в окне Гертцеля, например на середину. Те, которые обеспечивают один скалярный вывод, обычно не могут.

Вам также необходимо знать, где находится ваше окно Гертцеля относительно вашей временной оси.

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

Полное БПФ является избыточным, если вы знаете частоту вашего сигнала. Кроме того, Goertzel может быть настроен на частоту, непериодическую по длине БПФ, тогда как БПФ потребуется дополнительная интерполяция или заполнение нулями для непериодических в окне частот.

Комплексный Гертцель эквивалентен 1 бину ДПФ, который использует рекуррентность для векторов косинуса и синуса или коэффициентов твиндла БПФ.

hotpaw2
источник
Разве фазовая оценка нигде в пределах окна не имеет точно такой же точности, потому что вы просто добавили бы к фазовой оценке в начале окна, чтобы вычислить фазовую оценку в образце пределах окна ( быть началом окна)? к к = 0ωkkk=0
Олли Нимитало
Нет, потому что добавление wk приводит к другой фазе в конце окна, чем в начале для синусоиды с нецелым периодическим периодом в апертуре. Но 1-элементный ДПФ вычисляет одну круговую фазу в той же точке. Таким образом, все 3 значения будут разными. Но центральная фаза всегда связана с отношением нечетной / четной функции, независимо от того, что f0.
hotpaw2
Пытаюсь, но я не понимаю.
Олли Нимитало
Используйте косинус (нулевая фаза при k = 0), слегка подстройте частоту (на крошечное иррациональное число, но без изменения фазы при k = 0). DFT сообщает, что фаза изменилась! Попробуйте то же самое с косинусом точно по центру в k = N / 2. Без изменений при k = N / 2 для любого df. То же самое для греха или любой смеси. Центрирование контрольной точки фазы показывает меньшее изменение измеренной фазы с изменениями f0. например, ошибка частоты не способствует увеличению ошибок измерения фазы.
hotpaw2
1
Да, ошибка оценки фазы, которая меньше в центре окна, имеет смысл, если синусоида и фильтр Гоерцеля находятся на разных частотах. В этом случае оценка фазы, скажем, в конце окна смещена на постоянную величину, которая является произведением расстояния между центром и концом окна и разностью между частотами синусоидального фильтра и фильтра Гёртцела. Вычитание этого смещения дает ту же ошибку размера, что и для оценки центра, но требует знания частоты синусоиды.
Олли Нимитало
4

Если ваши сигналы не содержат шумов, вы можете определить пересечение нуля в обоих и определить частоту и относительную фазу.

Juancho
источник
3

Это зависит от того, что ваше определение «быстрый» является, насколько точным вы хотите , чтобы ваша оценка, хотите ли вы или фазу относительно ваших выборках, и как много шума есть на вашей функции и опорной синусоиды.ϕ

Один из способов сделать это - просто взять БПФ и просто посмотреть на корзину, ближайшую к . ωf(t)ω Однако это будет зависеть от того, что находится близко к центральной частоте бина.ω

Так:

  • Что вы подразумеваете под "быстро"?
  • Насколько точной вам нужна оценка?
  • Вы хотите (фаза относительно эталона) или фаза относительно начала отбора? Это имеет значение?ϕ
  • Каков уровень шума на каждом сигнале?

PS: я предполагаю, что вы имели в виду , а не .f(t)=Asin(ωt+ϕ)f(t)=Asin(ωx+ϕ)

Питер К.
источник
2

Начальная точка:
1) умножьте ваш сигнал и эталонную волну греха: = A⋅sin (ωt + ϕ) ⋅sin (ωt) = 0.5⋅A⋅ (cos (ϕ) - cos (2⋅ωt + ϕ) ) 2) найти интеграл по периоду : 3), который можно вычислить :
F(t)
T=π/ω
I(ϕ)=0TF(t)dt =0.5Acos(ϕ)T
ϕ
cos(ϕ)=I(t)/(0.5AT)

Задумайтесь:
как измерить А?
как определить в интервале ? (думать о «эталонный созе волны»)ϕ0..(2π)

Для дискретного сигнала измените интеграл на сумму и тщательно выберите T!

SergV
источник
1

Вы также можете сделать это (в нотации):

np.arctan( (signal*cos).sum() / (signal*sin).sum() ))

где сигнал - это сдвинутый по фазе сигнал, cos и sin - опорные сигналы, и вы генерируете аппроксимацию интеграла за определенное время путем суммирования по двум произведениям.

M529
источник
0

Это улучшение предложения @Kevin McGee об использовании одночастотного ДПФ с дробным индексом бина. Алгоритм Кевина не дает хороших результатов: хотя для половинок и целых бинов он очень точный, также близко к целым и половинкам, он также довольно хорош, но в противном случае ошибка может быть в пределах 5%, что, вероятно, неприемлемо для большинства задач. ,

Я предлагаю улучшить алгоритм Кевина, отрегулировав , то есть длину окна DFT, чтобы подходило как можно ближе к целому. Это работает, поскольку в отличие от FFT, DFT не требует, чтобы было степенью 2.NkN

Код ниже написан на Swift, но должен быть интуитивно понятен:

let f = 27.0 // frequency of the sinusoid we are going to generate
let S = 200.0 // sampling rate
let Nmax = 512 // max DFT window length
let twopi = 2 * Double.pi

// First, calculate k for Nmax, and then round it
var k = round(f * Double(Nmax) / S)

// The magic part: recalculate N to make k as close to whole as possible
// We also need to recalculate k once again due to rounding of N. This is important.
let N = Int(k * S / f)
k = f * Double(N) / S

// Generate the sinusoid
var r: [Double] = []
for i in 0..<N {
    let t = Double(i) / S
    r.append(sin(twopi * f * t))
}

// Compute single-frequency DFT
var R = 0.0, I = 0.0
let twopikn = twopi * k / Double(N)
for i in 0..<N {
    let x = Double(i) * twopikn
    R += r[i] * cos(x)
    I += r[i] * sin(x)
}
R /= Double(N)
I /= Double(N)

let amp = 2 * sqrt(R * R + I * I)
let phase = atan2(I, R) / twopi

print(String(format: "k = %.2f    R = %.8f    I = %.8f    A = %.8f    φ/2π = %.8f", k, R, I, amp, phase))
mojuba
источник
БПФ это просто способ эффективно рассчитать ДПФ. В современных библиотеках сила двух ограничений больше не существует. Если вам нужны только одно или два значения бина, лучше рассчитать их напрямую, как вы это сделали. Для одного чистого тона (действительного или сложного) необходимы только два значения бина, чтобы точно рассчитать частоту, фазу и амплитуду. См. Dsprelated.com/showarticle/1284.php . Математика довольно сложна, но есть ссылки на статьи, где объясняются производные. Линейная алгебра является предпосылкой для истинного понимания.
Седрон Доуг