Estou fazendo uma análise post-hoc de um conjunto de dados que consiste em uma série de quadros de filme que são contaminados por um artefato fortemente periódico. Gostaria de remover esse artefato dos meus quadros.
Para facilitar a plotagem, redimensionei minha matriz M
de valores de pixel [nframes, npixels]
e calculei a média de todos os valores de pixel para obter um vetor 1D m
. Veja como esse sinal se parece no domínio do tempo. Você pode ver a oscilação claramente na inserção ampliada.
Fiz então um periodograma tirando Fm = rfft(m)
e plotando abs(Fm)**2
contra frequência. Eu vejo um pico muito acentuado em ~ 1.5Hz:
Assim como a periodicidade temporal, também parece haver um componente espacial mais fraco para esse artefato, pois no valor exato da frequência de pico parece haver uma variação suave de fase no eixo x dos meus quadros, de modo que os pixels no direita tendem a ficar com pixels à esquerda:
Como abordagem de força bruta, tentei filtrar cada pixel no domínio do tempo usando um filtro de entalhe centrado em 1,5Hz. Usei um filtro Butterworth de ordem 4, com frequências críticas de 1,46 e 1,52Hz (não sou muito versado no design de filtros, por isso tenho certeza de que pode haver opções mais apropriadas).
Veja como é o sinal médio de pixel após a filtragem:
E o periodograma correspondente:
O filtro de entalhe faz um trabalho razoavelmente bom em reduzir o artefato, mas como ele basicamente se parece com um sinusóide estacionário puro, não posso deixar de pensar que poderia fazer melhor do que apenas atenuar essa parte do espaço de frequência.
Minha ideia inicial (muito ingênua) era fazer algo como:
- Obtenha a frequência, fase e amplitude da oscilação do espectro de Fourier para cada pixel do filme
- Reconstruir a oscilação no domínio do tempo
- Subtrair dos quadros de filme
Sei que isso não é algo que as pessoas costumam fazer, uma vez que a interferência geralmente não é tão espectralmente pura e temporalmente estacionária, mas me pergunto se isso pode fazer sentido no meu caso.
Dados
Pilha TIFF completa de 16 bits (~ 2 GB não compactada)
Versão de 8 bits dizimada espacialmente (~ 35 MB não compactada)
x
, então eu tomarFx = rfft(x)
, e obter o poder comoabs(Fx)**2
Respostas:
Sua solução proposta - calcular um senoide no domínio do tempo com base no pico da FFT e depois subtraí-lo - deve funcionar, mas há uma maneira mais fácil de fazer essencialmente a mesma coisa: modificar esse valor de pico na FFT e tomar o inverso transformar.
Portanto, para o seu vídeo rasterizado,
M[nframes, npixels]
você encontra o compartimento de frequência que contém o artefato e o achata sistematicamente (por exemplo, defina sua magnitude para a média de seus vizinhos) para cada pixel:Isso deve funcionar se o artefato tiver exatamente amplitude e frequência constantes, e sua frequência cair diretamente em um submúltiplo do comprimento da sequência (ou seja, os sinusóides representados pela FFT). Em geral, convém achatar uma ou duas caixas de cada lado
badbin
para lidar com um conjunto um pouco mais amplo de corrupções de banda estreita, por exemplo,Se você deseja restringir o componente removido de cada pixel a ter a mesma frequência e fase do artefato detectado na intensidade média, você pode remover apenas a projeção da
badbin
magnitude nessa fase, por exemplo,Observe que agora o componente resultante
badbin
sempre será 90 \ deg de fase (ortogonal) do globalbadbinphase
em cada pixel - qualquer componente de sinal exatamente nessa frequência e fase não pode ser separado do artefato.fonte