Transformada de Fourier. Os velozes e furiosos

Freqüentemente, ao desenvolver algoritmos, nos deparamos com o limite da complexidade computacional, que, ao que parece, é impossível superar. A transformação de Fourier tem complexidade , e uma versão rápida, proposta por volta de 1805 pela Casa 1 (e reinventada em 1965 por James Cooley e John Tukey) . Neste artigo, quero mostrar que você pode obter os resultados da conversão em tempo linear ou até alcançar dificuldade constante sob certas condições encontradas em problemas reais.


Quando me deparei com a tarefa de escrever um programa para analisar as funções de transferência de sistemas de som em tempo real, eu, como todo mundo, comecei pela conversão rápida. Tudo estava bem, mas com o grande tamanho da janela de tempo, a carga da CPU tornou-se indecentemente grande e algo teve que ser feito. Decidiu-se pausar e estudar a transformação novamente e, ao mesmo tempo, procurar maneiras de resolver o problema. Vamos voltar à transformação original de Joseph Fourier 2 :



Veremos com atenção o que está acontecendo aqui. Cada valor de saída no domínio da frequência é a soma de todos os valores de entrada do sinal multiplicado por . Para fazer cálculos, precisamos passar por todos os dados de entrada para cada valor de saída, ou seja, para cumprir aqueles operações.

Livre-se de n


Deixe-me lembrá-lo que, inicialmente, a tarefa era analisar dados de som em tempo real. Para fazer isso, a janela de tempo selecionada (essencialmente um buffer) de tamanho N é preenchida com dados com uma frequência f d correspondente à frequência de amostragem. No período T, os dados de entrada são convertidos da janela de tempo para a janela de frequência. Se você olhar para os números reais, N varia de 2 14 (16 384) a 2 16 (65 536) amostras (os valores são herdados da FFT, onde o tamanho da janela deve ser uma potência de dois). Tempo T = 80ms (12,5 vezes por segundo), que permite ver as alterações de maneira bastante conveniente e não sobrecarregar a CPU e a GPU. A frequência de amostragem f d é padrão e é de 48 kHz. Vamos calcular quantos dados na janela de tempo mudam entre dimensões. Durante o tempo T, ele entra no buffer amostras. Portanto, apenas 5% a 23% dos dados são atualizados na janela. Na pior das hipóteses, 95% (e no máximo 73%, o que também é muito!) Das amostras processadas cairão na conversão repetidas vezes, apesar de já terem sido processadas nas iterações anteriores.

O leitor atento naquele momento levantará a mão e dirá: “espere, mas e o coeficiente ? Afinal, a cada nova transformação, os mesmos dados serão localizados nas novas posições da série e, como resultado, terão coeficientes diferentes? ” Para cada cinco cuidados, lembremos de um detalhe importante na transformação que geralmente é esquecida. No estudo dos valores das funções no intervalo de 0 a t, a função é considerada periódica, o que permite que você mude sem esforço a função para a esquerda ou direita no tempo. Como resultado, temos todo o direito de não inserir um novo valor no final e excluir o valor antigo desde o início, mas substituir ciclicamente os dados no buffer.

Para maior clareza, você pode escrever em forma de tabela como o buffer será alterado:
t = 0f (0)f (1)f (2)f (3)f (4)f (5)f (6)f (7)f (8)f (9)
t = 1f (10)f (1)f (2)f (3)f (4)f (5)f (6)f (7)f (8)f (9)
t = 2f (10)f (11)f (2)f (3)f (4)f (5)f (6)f (7)f (8)f (9)
t = 3f (10)f (11)f (12)f (3)f (4)f (5)f (6)f (7)f (8)f (9)
t = 4f (10)f (11)f (12)f (13)f (4)f (5)f (6)f (7)f (8)f (9)

Você pode escrever como a transformação no tempo muda de t 1 para t 2 :


Valor é o resultado da conversão anterior e a complexidade do cálculo não depende do tamanho da janela de tempo e, portanto, é constante. Como resultado, a complexidade da conversão será * porque tudo o que resta para nós é atravessar a janela de frequência uma vez e aplicar as alterações nas amostras T que mudaram ao longo do tempo. Gostaria também de chamar sua atenção para o fato de que as probabilidades pode ser calculado antecipadamente, o que fornece um ganho adicional de produtividade, e apenas duas operações permanecem no ciclo: subtraindo números reais e multiplicando um número real por um complexo, na prática, essas duas operações são simples e baratas.

Para completar a imagem, resta apenas indicar o estado inicial, mas aqui tudo é simples:


* - é claro, a complexidade final de toda a transformação permanecerá tão , mas será executado gradualmente, em n iterações, enquanto o buffer for atualizado. - essa é a complexidade da atualização dos dados, mas é exatamente isso que precisamos (ao usar a FFT, a complexidade de cada transformação )

Mas e se você se aprofundar. Ou se livrar do segundo n


Quero fazer uma reserva imediata de que as próximas etapas serão aplicáveis ​​apenas se você não planeja executar a transformação inversa para o resultado (para corrigir o sinal ou obter uma resposta de impulso). Para começar, quero lembrá-lo de que, como resultado da conversão, obtemos uma matriz de dados simétrica, o que imediatamente nos permite reduzir pela metade o número de conversões.

Agora vamos analisar o conjunto de dados resultante, dadas as condições do problema. Temos um conjunto de números complexos, cada um dos quais descreve a amplitude e a fase das oscilações em uma frequência específica. A frequência pode ser determinada pela fórmula: para . Vamos avaliar a etapa da janela de frequência em nossos dados: Para N = 2 14 : 2,93 Hz (e para 2 16 : 0,73 Hz). Assim, na faixa de 1 kHz a 2 kHz, obtemos 341 resultados. Tente avaliar de forma independente quantos dados estarão na faixa de 8 kHz a 16 kHz para N = 65536. Muito, certo? Muito! Precisamos de tantos dados? Obviamente, nos problemas de exibição das características de frequência dos sistemas de som, a resposta é não. Por outro lado, para análise na região de baixa frequência, um pequeno passo é muito útil. Não se esqueça que ainda há um cronograma pela frente que deve esses volumes ( ) converta para um formato legível por humanos (média, spline ou suavização) e exiba-os na tela. E em altas frequências, mesmo tendo uma tela 4K e exibindo o gráfico no modo de tela cheia com o eixo de frequência logarítmica, o tamanho da etapa rapidamente se tornará muito menor que 1 pixel.

Pela experiência, você pode descobrir que basta ter apenas 48 pontos por oitava e ter os dados um pouco mais suaves e com média mais alta, sugiro parar em 96. Na faixa de frequência de áudio de 20 Hz a 20 kHz, é fácil contar apenas 10 oitavas: 20, 40, 80 , 160, 320, 640, 1280, 2560, 5120, 10240, 20480, cada um dos quais pode ser dividido em um determinado número de sub-bandas (não esqueça que a partição deve ser feita geometricamente e não aritmeticamente); portanto, é mais do que suficiente realizar a conversão apenas para 960 frequências para obter Performan que em 16 ... 65 vezes menor do que a versão original.

Assim, combinando as duas abordagens, obtemos a complexidade constante do algoritmo de atualização de dados .

Mel ao quadrado e uma colher de alcatrão


Agora podemos dizer com segurança que pela complexidade chegamos à complexidade usando dois hacks simples da vida:

  • Depois de analisar o problema, percebemos que os dados são adicionados gradualmente, e o período de atualização completa da janela de tempo é muito superior ao período de transformações e calculamos a diferença da transformação de Fourier.
  • passou da etapa aritmética na janela de frequência para limitado apenas pelos valores especificados, o que pode reduzir drasticamente o número de conversões.

Mas, é claro, a vida seria realmente um conto de fadas, se não um. A aplicação dessas duas abordagens nos permitiu realmente descarregar a CPU, de modo que adivinhar que ela calcula a transformação de Fourier e exibe os resultados na tela, mesmo com foi dificil Mas a punição não ficou esperando quando seus sinais na realidade não são periódicos (e isso é necessário para obter os resultados de conversão corretos) e não é possível selecionar o tamanho de janela apropriado, torna-se necessário o uso de várias funções da janela, o que não permite mais que você use a primeira etapa completamente. A prática demonstrou que o uso das funções da janela é crítico no estudo de sinais com frequência menor que . Em altas frequências, o número de períodos que caem na janela de tempo atenua significativamente as distorções decorrentes da presença de um intervalo de primeira ordem (entre f (0) ef (N-1)) na função original.

No final, eu também recusei o segundo passo e voltei para a FFT, porque o ganho nessa tarefa já era pequeno.

Em conclusão


A primeira abordagem pode ser aplicada se seus dados forem de natureza periódica pronunciada e precisarem ser analisados ​​ao longo do tempo usando uma grande janela de tempo, que, lembro-me, não precisa ser do grau 2, ou seja, qualquer número natural.
A segunda abordagem é aplicável (mesmo levando em consideração as funções da janela) se apenas um pequeno conjunto de frequências for analisado nos dados.

Infelizmente, para mim neste problema, isso permaneceu apenas um pouco de entretenimento matemático, mas espero que o inspire a estudar outros algoritmos nos feriados em termos de alterações nos dados de entrada ao longo do tempo :)

Literatura


  1. Cooley - algoritmo Tukey FFT
  2. E.A. Vlasova Rows. Editora da MSTU com o nome de N.E. Bauman. Moscovo 2002

Imagem retirada do mangá de Michio Shibuya. “EMOCIONANTE MATEMÁTICA. Análise de Fourier

Source: https://habr.com/ru/post/pt434734/


All Articles