Como posso calcular um espectro de energia com espaçamento de log?

Eu gostaria de calcular um espectro de potência no qual as frequências são espaçadas logaritmicamente.

No método de Welch, há uma troca entre a resolução de frequência do espectro de potência resultante e o número de médias (isto é, erro no resultado). Eu gostaria que esse trade-off fosse dinâmico, ou seja, faça menos médias para pontos de baixa frequência para ter uma resolução mais precisa em baixa frequência.

Existe uma maneira padrão de fazer isso?

Suponho que uma maneira seria inicialmente fazer pwelchcom uma resolução muito alta (baixo número de médias) e, em seguida, refinar o espectro resultante usando o binar logarítmico.

Encontrei um artigo que aborda essa questão diretamente:

• Michael Tröbs e Gerhard Heinzel, " Melhoria da estimativa de espectro a partir de séries temporais digitalizadas em um eixo de frequência logarítmica ", em Measurement , vol 39 (2006), pp 120-129. (Ou pré-impressão gratuita .)

As primeiras figuras do artigo ilustram bem o problema que esse algoritmo resolve, e as referências contêm uma bibliografia útil de outras abordagens (transformação Q constante, transformação de Fourier temperada, artigo de pesquisa etc.).

Sua abordagem não é refazer a saída de uma estimativa do espectro de potência baseada em FFT, mas calcular apenas a transformada de Fourier discreta nas frequências de interesse (espaçadas logaritmicamente). Para cada frequência a ser estimada, eles basicamente implementam o algoritmo de Welch, mas com um comprimento de transformação (e, portanto, também, número de médias) escolhido especificamente para cada frequência. O cálculo de cada compartimento de frequência usa toda a série temporal, mas segmentada de maneira diferente. Os resultados têm a propriedade desejável de que a resolução (largura do compartimento) seja uma função suave da frequência e os resultados podem ser calibrados como uma densidade espectral de potência ou um espectro de potência.

Implementação do Matlab aqui: https://github.com/tobin/lpsd

Divulgação: Os autores deste artigo estão na mesma instituição que eu.

Nesse caso, eu usaria um método de mínimos quadrados para calcular a frequência de alguma lista de valores conhecida. O método mais comum é o método Lomb. Funciona de maneira bastante semelhante a uma FFT ou DFT, mas calcula apenas a frequência em frequências determinadas e pode lidar com dados ausentes, caso isso seja um problema. A ideia é a seguinte:

1. $w$
2. $w$$t_j$$X_j$

$P_x(\omega) = \frac{1}{2} \left( \frac { \left[ \sum_j X_j \cos \omega ( t_j - \tau ) \right] ^ 2} { \sum_j \cos^2 \omega ( t_j - \tau ) } + \frac {\left[ \sum_j X_j \sin \omega ( t_j - \tau ) \right] ^ 2} { \sum_j \sin^2 \omega ( t_j - \tau ) } \right)$

Observe que isso não será tão bem quanto um FFT, portanto, eu o faria apenas se o número de frequências desejadas for muito menor do que o FFT, que seria necessário para coletar todos os dados.

Caso contrário, poderia-se fazer um método de interpolação ou qualquer outra re-amostragem de uma FFT ou DFT.