Existe um algoritmo para calcular a fase para uma única frequência?

Se você tem uma função $f(t)=A \cdot \sin(\omega t+\phi)$ e faz referência à onda $\sin(\omega x)$ qual seria um algoritmo rápido para calcular $\phi$ ?

Eu estava olhando para o algoritmo de Goertzel , mas ele não parece lidar com a fase?

Use um DFT na frequência específica. Em seguida, calcule a amplitude e a fase das partes real / imag. Ele fornece a fase referida ao início do tempo de amostragem.

Em uma FFT 'normal' (ou uma DFT calculada para todos os N harmônicos), você normalmente calcula a frequência com f = k * (sample_rate) / N, onde k é um número inteiro. Embora possa parecer sacrílego (especialmente para membros da Igreja do Inteiro Inteiro), você pode realmente usar valores não inteiros de k ao fazer uma única DFT.

Por exemplo, suponha que você tenha gerado (ou obtido) N = 256 pontos de uma onda senoidal de 27 Hz. (digamos, sample_rate = 200). Suas frequências 'normais' para uma FFT de 256 pontos (ou D N de ponto N) corresponderiam a: f = k * (sample_rate) / N = k * (200) / 256, em que k é um número inteiro. Mas um 'k' não inteiro de 34,56 corresponderia a uma frequência de 27 Hz., Usando os parâmetros listados acima. É como criar um 'compartimento' de DFT que esteja exatamente centralizado na frequência de interesse (27 Hz.). Alguns códigos C ++ (compilador DevC ++) podem ter a seguinte aparência:

#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: Espero que o texto acima traduza bem o stackoverflow - alguns deles podem ser contornados)

O resultado do exposto acima é uma fase de -twopi / 4, conforme mostrado nos pontos reais gerados (e o amplificador é duplicado para refletir a frequência pos / neg).

Algumas coisas a serem observadas - eu uso o cosseno para gerar a forma de onda de teste e interpretar os resultados - você precisa ter cuidado com isso - a fase é referenciada no tempo = 0, que é quando você começou a amostrar (ou seja, quando você coletou r [0] ) e cosseno é a interpretação correta).

O código acima não é elegante nem eficiente (por exemplo: use tabelas de consulta para os valores sin / cos, etc.).

Seus resultados ficarão mais precisos quando você usar N maior, e há um pouco de erro devido ao fato de a taxa de amostragem e N acima não serem múltiplos um do outro.

Obviamente, se você quiser alterar sua taxa de amostragem, N ou f, precisará alterar o código e o valor de k. Você pode colocar uma bandeja DFT em qualquer lugar da linha de frequência contínua - apenas certifique-se de usar um valor de k que corresponde à frequência de interesse.

O problema pode ser formulado como problema dos mínimos quadrados (não linear):

$F(\phi)$$\phi$

A derivada é muito simples:

$F''(\phi)$

Obviamente, a função objetivo acima tem vários mínimos por causa da periodicidade, portanto, alguns termos de penalidade podem ser adicionados para discriminar outros mínimos (por exemplo, adicionando à equação do modelo). Mas acho que a otimização só vai convergir para os mínimos mais próximo e você pode atualizar o resultado subtraindo . 2 π k , k ∈ $\phi^{2}$$2\pi k, k\in N$

Existem várias formulações diferentes do algoritmo de Goertzel. As que fornecem 2 variáveis ​​de estado (ortogonais ou próximas a) ou uma variável de estado complexa, como saídas possíveis, muitas vezes podem ser usadas para calcular ou estimar a fase com referência a algum ponto da janela de Goertzel, como o meio. Aqueles que fornecem apenas uma saída escalar normalmente não podem.

Você também precisará saber onde está sua janela do Goertzel em relação ao seu eixo temporal.

Se o seu sinal não for exatamente inteiro periódico na janela do Goertzel, a estimativa da fase em torno de um ponto de referência no meio da janela poderá ser mais precisa do que a fase de referência no início ou no final.

Uma FFT completa é um exagero se você souber a frequência do seu sinal. Além disso, um Goertzel pode ser sintonizado para uma frequência não periódica no comprimento da FFT, enquanto uma FFT precisará de interpolação adicional ou preenchimento zero para frequências não periódicas na janela.

Um Goertzel complexo é equivalente a 1 bin de uma DFT que usa uma recorrência para os vetores de base cosseno e senoidal ou fatores de variação da FFT.

Se seus sinais estiverem livres de ruído, você poderá identificar cruzamentos zero em ambos e determinar a frequência e a fase relativa.

Isso depende de qual é a sua definição de "rápido", quão precisa você deseja sua estimativa, se deseja ou a fase relativa às suas amostragens e quanto ruído há em sua função e onda senoidal de referência.$\phi$

Uma maneira de fazer isso é pegar a FFT de e apenas olhar para a lixeira mais próxima de . $f(t)$$\omega$ No entanto, isso dependerá de estar próximo da frequência do centro da bandeja.$\omega$

Assim:

• O que você quer dizer com "rápido"?
• Qual é a precisão da estimativa?
• $\phi$
• Qual é o nível de ruído em cada sinal?

PS: Eu estou assumindo que você quis dizer , em vez de .$f(t)=A\sin(ωt+ϕ)$$f({\Huge t})=A\sin(ω{\Huge x}+ϕ)$

Ponto de partida:
1) multiplique seu sinal e a onda de pecado de referência: = A⋅sin (ωt + ϕ) ⋅sin (ωt) = 0.5⋅A⋅ (cos (ϕ) - cos (2⋅ωt + ϕ) ) 2) encontre integral no período : 3) você pode calcular :
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)$

Pense em:
como medir A?
como determinar no intervalo ? (pense em " onda cos de referência ")0 .. ( 2 ⋅ ¸ $\phi$$0..(2 \cdot \pi)$

Para sinal discreto, altere a integral para somar e escolha cuidadosamente T!

Você também pode fazer isso (em notação numpy):

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

onde o sinal é o seu sinal de mudança de fase, cos e sin são os sinais de referência, e você gera uma aproximação de uma integral ao longo de um certo tempo através da soma dos dois produtos.

Essa é uma melhoria na sugestão de @Kevin McGee de usar um DFT de frequência única com um índice de escaneamento fracionário. O algoritmo de Kevin não produz ótimos resultados: embora com meias caixas e caixas inteiras seja muito preciso, também próximo aos totais e metades também é muito bom, mas, caso contrário, o erro pode estar dentro de 5%, o que provavelmente não é aceitável para a maioria das tarefas .

Sugiro melhorar o algoritmo de Kevin ajustando , ou seja, o comprimento da janela DFT para que fique o mais próximo possível de um todo. Isso funciona, pois, diferentemente da FFT, a DFT não exige que seja uma potência de 2.$N$$k$$N$

O código abaixo está no Swift, mas deve ser intuitivamente claro:

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