Como posso identificar e remover automaticamente uma frequência fundamental e todos os seus harmônicos?

Eu tenho alguns dados de microscopia que estão contaminados por um artefato de batimento cardíaco que eu gostaria de remover. Os dados consistem em uma grande série temporal de imagens capturadas em ~ 60Hz.

Aqui está um pequeno exemplo de clipe em formato GIF:

Peguei a intensidade média de pixels ao longo do tempo e calculei o periodograma usando o método de Welch:

Como você pode ver, há um pico acentuado em ~ 1,8 Hz, que provavelmente corresponde à freqüência cardíaca (~ 108 batimentos / min). Há também vários picos harmônicos em múltiplos inteiros de 1,8Hz. É provável que a freqüência cardíaca exata varie de conjunto de dados para conjunto de dados, mas posso especificar um intervalo biologicamente plausível, conforme mostrado pela área sombreada no periodograma.

O que eu gostaria de poder fazer é:

1. Detecta automaticamente a frequência fundamental correspondente ao batimento cardíaco e todos os seus harmônicos
2. Filtre os dados para remover os harmônicos fundamentais e todos.

No momento, posso resolver o ponto 1 de maneira muito grosseira, encontrando o maior pico no periodograma e depois multiplicando-o por $1, 2, ..., N$ Onde $N$ é o número estimado de picos harmônicos, mas tenho certeza de que deve haver um método melhor que esse hack.

Em relação ao ponto 2, deparei-me com essa pergunta, que menciona o uso de um filtro de pente para remover um fundamental e todos os seus harmônicos. Esse é o melhor método para usar? Uma consideração importante é que terei de aplicar o filtro a cada série temporal de pixels em uma grande matriz, portanto, um método computacionalmente eficiente seria altamente desejável.

Dados de exemplo

• Fiz upload de uma série de tempo 1D com exemplos de valores médios de pixel em .csvformato aqui .
• Também adicionei uma amostra de 1000 quadros (dizimados espacialmente 2: 1) em um .matarquivo compactado aqui

Seu método não é tão ruim para uma primeira tentativa.

No entanto, o método a seguir tende a funcionar melhor:

1. Procure maxima local
2. Verifique os máximos próximos (espaçamento de 2 ou 3 posições) e mescle-os
3. Crie algumas hipóteses sobre a frequência fundamental. Atualmente, você assume que o pico mais alto é a frequência fundamental, que é uma hipótese. Você também deve verificar a possibilidade de que o pico mais alto seja o primeiro harmônico, ou seja, que exista um pico menor na metade da frequência do pico principal. Você também pode ter outros casos a considerar, usando seu conhecimento do problema em questão (interferência? Batimento cardíaco irregular?).
4. Assumindo cada uma dessas hipóteses, encontre a frequência do solo ajustando a parábola a cada pico harmônico. Cada pico produzirá uma estimativa ligeiramente diferente devido ao ruído, mas esses erros não estão correlacionados e estão fora da média. Uma das hipóteses levará a ajustes muito melhores, escolha a frequência do solo prevista por essa hipótese.
5. Usando a frequência do solo que você encontrou na etapa 4, especifique a parábola ao redor de cada pico para estimar a altura do pico. Observe que os picos provavelmente caem entre duas caixas.
6. Agora você tem a localização e a força do fundamental e de seus harmônicos, mas não da fase. Provavelmente é mais fácil encontrar a fase do fundamental, subtrair isso, encontrar a fase do primeiro harmônico etc.

A principal razão pela qual isso funciona melhor é a etapa 4. Qualquer hipótese errada para uma frequência de terra falhará terrivelmente, enquanto você tenta ajustar picos ao redor da localização das harmônicas previstas. Digamos que você tenha um pico em 2 Hz. Essa pode ser a frequência do solo ou o primeiro harmônico. Quando você testa a "hipótese da primeira harmônica", ou seja, se a frequência do solo é realmente 1 Hz, você ajusta as parábolas aos dados em torno de 1,2,3,4,5 ... Hz. Se essa hipótese estiver errada, você recebe lixo em torno de 1,3,5 Hz. Se estiver certo, você poderá encontrar picos próximos a 1,1 Hz, 2,2, 3,3, 4,4 e 5,6 - o que sugere que a frequência real do solo é de 1,12 Hz.

Você está procurando subtração espectral iterativa . Aqui estão algumas informações gerais de um livro de Análise de conteúdo de Alexander Lerch.

Eu recomendaria uma correlação automática no seu periodograma. Você pode construir filtros de pico de entalhe ou negativos com base em múltiplos da frequência fundamental que a correlação produz.

Esse código me ajudou a criar um bom gráfico de correlação automática (cepstrum é outro bom método quando os harmônicos têm mais poder do que o fundamental) http://note.sonots.com/SciSoftware/Pitch.html

1. Automatically detect the fundamental frequency corresponding to the heartbeat, and all of its harmonics

Você pode fazer uma média local na DFT para encontrar os batimentos cardíacos. Se um determinado ponto desse grupo for maior que thresholde o maior valor em alguma faixa em torno dele, será um batimento cardíaco ou harmônico.

2. Filter the data so as to remove the fundamental and all harmonics.

Você pode simplesmente usar os índices de batimentos cardíacos localizados anteriormente e costurá-los linearmente. Se isso não for muito grosseiro para você.

Qual é o objetivo desta parte da operação? Esta é uma análise médica que requer a ausência do batimento cardíaco para procurar outros dados (talvez) ocultos? Ou isso é apenas para estética?

Eu imagino que um filtro de pente será um pouco incontrolável. Isso exigirá um ajuste cuidadoso do parâmetro de feedback. Pode ser viável.

EDIT: Você precisa de um algoritmo de detecção de afinação? Escrevi um desses há alguns anos ignorando os termos (próximos a) do DC e encontrando o primeiro máximo local. Depois, usando os compartimentos de ambos os lados, eu poderia usar uma interpolação quadrática e encontrar a localização da frequência máxima com mais precisão do que a resolução permitida pelas frequências dos próprios compartimentos.

Seria mais preciso utilizar harmônicos? Provavelmente isso apenas adicionaria uma oportunidade para erros. Embora você possa tentar obter um valor para o fundamental com esse método e encontrar o máximo em dobro, triplo, etc, e depois usar um método semelhante ao anterior:

1. encontre o máximo local 2. interpolado quadrático para encontrar o valor entre caixas para frequência.

Se você vetar a etapa quadrática e apenas der o máximo local, não obterá uma frequência precisa e subir os harmônicos ajudaria.