Removendo um Artefato Senoidal de um Conjunto de Quadros de Filme

7

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 Mde 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.

insira a descrição da imagem aqui

Fiz então um periodograma tirando Fm = rfft(m)e plotando abs(Fm)**2contra frequência. Eu vejo um pico muito acentuado em ~ 1.5Hz:

insira a descrição da imagem aqui

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:

insira a descrição da imagem aqui

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: insira a descrição da imagem aqui

E o periodograma correspondente: insira a descrição da imagem aqui

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:

  1. Obtenha a frequência, fase e amplitude da oscilação do espectro de Fourier para cada pixel do filme
  2. Reconstruir a oscilação no domínio do tempo
  3. 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)

ali_m
fonte
A partir de uma série de quadros de filme, você pode esclarecer com mais clareza como está gerando exatamente o PSD?
Tarin Ziyaee
@ user4619 de forma muito grosseira - para cada quadro Acabei calculado o valor médio do pixel para gerar um vetor x, então eu tomar Fx = rfft(x), e obter o poder comoabs(Fx)**2
ali_m
Você tem um quadro 2-D e gera um vetor 1-D médio. Junto x? Ao longo de y?
Tarin Ziyaee
@ user4619 junto ambos x e y - Eu remodelar o meu filme em um nframes por variedade npixels, em seguida, média de todas pixels
ali_m
Ok, obrigado por esse detalhe - é importante na análise. Adicione esta informação à sua publicação.
Tarin Ziyaee

Respostas:

1

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:

import numpy as np
nframes, npixels = np.shape(M)
# Identify the bin containing the sinusoidal artifact
# Use the average intensity for each image
m = np.mean(M, axis=1)
# Calculate the FFT
Fm = np.fft.rfft(m)
# Find the largest bin away from the low-frequency region
lowfreq = 100  # or something
badbin = lowfreq + np.argmax(Fm[lowfreq:]**2)

# Now adjust the amplitude of that bin in the FFT of each pixel
for pixel in range(npixels):
   Fpix = np.fft.rfft(M[:, pixel])
   # Scale magnitude of artifact bin to be the mean of its neighbors
   Fpix[badbin] *= np.mean(np.absolute(Fpix[[badbin-1, badbin+1]]))/np.absolute(Fpix[badbin])
   # Rewrite the time sequence of that pixel
   M[:,pixel] = np.fft.irfft(Fpix)

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 badbinpara lidar com um conjunto um pouco mais amplo de corrupções de banda estreita, por exemplo,

# ...

   # Scale magnitude of artifact binS to be the mean of neighbors
   spread = 3  # flatten bins from (badbin - (spread-1)) to (badbin + (spread-1))
   # target value for new bins
   targetmag = np.mean(np.absolute(Fpix[[badbin-spread, badbin+spread]]))
   bins = range(badbin - (spread-1), badbin + spread)
   Fpix[bins] *= targetmag/np.abs(Fpix[bins])
   # ...

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 badbinmagnitude nessa fase, por exemplo,

badbinphase = np.angle(Fm[badbin])
# ...

   Ncomponent = np.abs(Fpix[badbin])*np.cos(np.angle(Fpix[badbin]) - badbinphase)
   Fpix[badbin] -= Ncomponent * np.exp(0+1j * badbinphase)
   # ...

Observe que agora o componente resultante badbinsempre será 90 \ deg de fase (ortogonal) do global badbinphaseem cada pixel - qualquer componente de sinal exatamente nessa frequência e fase não pode ser separado do artefato.

dpwe
fonte
Não é essencialmente isso que já estou fazendo com o filtro de entalhe? Ainda acho que não é o ideal, já que essa abordagem não leva em conta as informações sobre a fase da oscilação que estou tentando remover, que é constante ao longo do tempo. Parece-me que deve ser possível remover seletivamente o artefato sem afetar o sinal 'genuíno' que cai na mesma faixa de frequência.
ali_m
Não foi exatamente o que você fez: primeiro, é um entalhe muito mais estreito; segundo, ao contrário do filtro Butterworth, ele não causa distorção de fase. Se você subtrair as FFTs antes e depois da modificação, obterá um único componente diferente de zero com alguma amplitude e a fase do pico original. Este é o sinal que estamos removendo no domínio do tempo, ou seja, um sinusóide de amplitude constante exatamente na frequência e fase do pico espectral. Se o sinal subjacente tiver energia nessa região, ele aparecerá como ruído na estimativa, mas deverá desaparecer.
dpwe
Sei que você pode ter significado impor a mesma fase em todos os pixels, então acrescentei isso à minha resposta. Porém, não é mágico - você não pode separar o componente em fase do sinal e do ruído; portanto, você sempre fica com 90 \ deg residual do artefato estimado.
dpwe