Existe-t-il un algorithme pour calculer la phase pour une seule fréquence?

Si vous avez une fonction $f(t)=A \cdot \sin(\omega t+\phi)$ , et l' onde sin de référence $\sin(\omega x)$ ce serait un algorithme rapide pour calculer $\phi$ ?

Je regardais l' algorithme de Goertzel , mais il ne semble pas traiter de phase?

Utilisez un DFT à la fréquence spécifique. Calculez ensuite l'amplitude et la phase à partir des parties réelles / imagées. Il vous donne la phase référencée au début du temps d'échantillonnage.

Dans une FFT 'normale' (ou une DFT calculée pour toutes les N harmoniques), vous calculez généralement la fréquence avec f = k * (sample_rate) / N, où k est un entier. Bien que cela puisse sembler sacrilège (en particulier pour les membres de l'Église du Entier Entier), vous pouvez réellement utiliser des valeurs non entières de k lorsque vous effectuez une seule DFT.

Par exemple, supposons que vous ayez généré (ou obtenu) N = 256 points d'une onde sinusoïdale de 27 Hz. (disons, sample_rate = 200). Vos fréquences «normales» pour une FFT à 256 points (ou DFT à N points) correspondraient à: f = k * (taux_échantillon) / N = k * (200) / 256, où k est un entier. Mais un «k» non entier de 34,56 correspondrait à une fréquence de 27 Hz., En utilisant les paramètres énumérés ci-dessus. C'est comme créer un «bac» DFT qui est exactement centré sur la fréquence d'intérêt (27 Hz.). Certains codes C ++ (compilateur DevC ++) peuvent ressembler à ceci:

#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: j'espère que ce qui précède se traduit bien par stackoverflow - une partie pourrait se terminer)

Le résultat de ce qui précède est une phase de -twopi / 4, comme indiqué dans les points réels générés (et l'ampli est doublé pour refléter la fréquence pos / neg).

Quelques points à noter - j'utilise le cosinus pour générer la forme d'onde de test et interpréter les résultats - vous devez être prudent à ce sujet - la phase est référencée à time = 0, c'est-à-dire lorsque vous avez commencé l'échantillonnage (c'est-à-dire lorsque vous avez collecté r [0] ), et cosinus est la bonne interprétation).

Le code ci-dessus n'est ni élégant ni efficace (par exemple: utilisez des tables de recherche pour les valeurs sin / cos, etc.).

Vos résultats seront plus précis à mesure que vous utiliserez un N plus grand, et il y a un peu d'erreur en raison du fait que la fréquence d'échantillonnage et N ci-dessus ne sont pas des multiples l'un de l'autre.

Bien sûr, si vous souhaitez modifier votre fréquence d'échantillonnage, N ou f, vous devrez modifier le code et la valeur de k. Vous pouvez plonger dans un bac DFT n'importe où sur la ligne de fréquence continue - assurez-vous simplement que vous utilisez une valeur de k qui correspond à la fréquence d'intérêt.

Le problème peut être formulé comme un problème de moindres carrés (non linéaire):

$F(\phi)$$\phi$

Le dérivé est très simple:

$F''(\phi)$

De toute évidence, la fonction objectif ci-dessus a plusieurs minima en raison de la périodicité, donc un terme de pénalité peut être ajouté pour discriminer d'autres minima (par exemple, en ajoutant à l'équation du modèle). Mais je pense que l'optimisation va simplement converger vers les minima le plus proche et vous pouvez mettre à jour le résultat soustrayant . 2 π k , k ∈ $\phi^{2}$$2\pi k, k\in N$

Il existe plusieurs formulations différentes de l'algorithme de Goertzel. Celles qui fournissent 2 variables d'état (orthogonales ou proches de), ou une variable d'état complexe, car les sorties possibles peuvent souvent être utilisées pour calculer ou estimer la phase en référence à un certain point de la fenêtre de Goertzel, comme le milieu. Ceux qui fournissent une seule sortie scalaire ne le peuvent généralement pas.

Vous devrez également savoir où se trouve votre fenêtre Goertzel par rapport à votre axe temporel.

Si votre signal n'est pas exactement un entier périodique dans votre fenêtre Goertzel, l'estimation de phase autour d'un point de référence au milieu de la fenêtre peut être plus précise que la phase de référence au début ou à la fin.

Une FFT complète est exagérée si vous connaissez la fréquence de votre signal. De plus, un Goertzel peut être réglé sur une fréquence non périodique dans la longueur de la FFT, tandis qu'une FFT aura besoin d'une interpolation supplémentaire ou d'un remplissage nul pour les fréquences non périodiques dans la fenêtre.

Un Goertzel complexe équivaut à 1 casier d'un DFT qui utilise une récurrence pour les vecteurs de base cosinus et sinus ou les facteurs de torsion FFT.

Si vos signaux sont exempts de bruit, vous pouvez identifier les passages par zéro dans les deux et déterminer la fréquence et la phase relative.

Cela dépend de votre définition de "rapide", de la précision de votre estimation, de la valeur de ou de la phase par rapport à vos échantillonnages, et du niveau de bruit sur votre fonction et votre onde sinusoïdale de référence.$\phi$

Une façon de le faire est de simplement prendre la FFT de et de regarder simplement le bac le plus proche de . $f(t)$$\omega$ Cependant, cela dépendra du fait que est proche de la fréquence centrale du bac.$\omega$

Donc:

• Qu'entendez-vous par «rapide»?
• Dans quelle mesure avez-vous besoin de l'estimation?
• Voulez-vous (phase relative à la référence) ou phase relative au début de l'échantillonnage? Est-ce que ça importe?$\phi$
• Quel est le niveau de bruit sur chaque signal?

PS: Je suppose que vous vouliez dire , plutôt que .$f(t)=A\sin(ωt+ϕ)$$f({\Huge t})=A\sin(ω{\Huge x}+ϕ)$

Point de départ:
1) multipliez votre signal et l'onde sinusoïdale de référence: = A⋅sin (ωt + ϕ) ⋅sin (ωt) = 0,5⋅A⋅ (cos (ϕ) - cos (2⋅ωt + ϕ) ) 2) trouver l'intégrale sur la période : 3) vous pouvez calculer :
T = π / ω I ( ϕ ) = ∫ T 0 F ( t ) d t = 0,5 ⋅ A ⋅ c o s ( ϕ ) ⋅ T ϕ c o s ( ϕ ) = I ( t ) / ( 0,5 ⋅ A ⋅ T $F(t)$
$T= \pi /\omega$
$I(\phi) = \int_0^T F(t)dt\ = 0.5⋅A⋅cos(ϕ) \cdot T$
$\phi$
$cos(\phi) = I(t)/(0.5 \cdot A \cdot T)$

Pensez à:
comment mesurer A?
comment déterminer dans l' intervalle ? (pensez à " onde cos de référence ")0 .. ( 2 ⋅ π $\phi$$0..(2 \cdot \pi)$

Pour un signal discret, changez l'intégrale pour additionner et choisissez soigneusement T!

Vous pouvez également faire cela (en notation numpy):

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

où signal est votre signal déphasé, cos et sin sont les signaux de référence, et vous générez une approximation d'une intégrale sur un certain temps via la sommation sur les deux produits.

Il s'agit d'une amélioration par rapport à la suggestion de @Kevin McGee d'utiliser un DFT à fréquence unique avec un indice de groupe fractionnaire. L'algorithme de Kevin ne donne pas d'excellents résultats: alors qu'il est très précis dans les demi-bacs et les bacs entiers, il est également assez proche des entiers et des moitiés, mais sinon l'erreur peut être inférieure à 5%, ce qui n'est probablement pas acceptable pour la plupart des tâches .

Je suggère d'améliorer l'algorithme de Kevin en ajustant , c'est-à-dire la longueur de la fenêtre DFT afin que se rapproche le plus possible d'un tout. Cela fonctionne car contrairement à la FFT, la DFT ne nécessite pas que soit une puissance de 2.$N$$k$$N$

Le code ci-dessous est en Swift, mais devrait être intuitivement clair:

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))