¿Cómo puedo identificar y eliminar automáticamente una frecuencia fundamental y todos sus armónicos?

Tengo algunos datos de microscopía que están contaminados por un artefacto de latido que me gustaría eliminar. Los datos consisten en una gran serie temporal de imágenes capturadas a ~ 60Hz.

Aquí hay un pequeño clip de ejemplo en forma GIF:

Tomé la intensidad de píxel promedio a lo largo del tiempo y calculé el periodograma usando el método de Welch:

Como puede ver, hay un pico agudo a ~ 1.8Hz que probablemente corresponderá a la frecuencia cardíaca (~ 108 latidos / min). También hay un montón de picos armónicos en múltiplos enteros de 1.8Hz. Es probable que la frecuencia cardíaca exacta varíe de un conjunto de datos a otro, pero puedo especificar un rango biológicamente plausible como se muestra en el área sombreada en el periodograma.

Lo que me gustaría poder hacer es:

1. Detecta automáticamente la frecuencia fundamental correspondiente al latido del corazón y todos sus armónicos.
2. Filtre los datos para eliminar el fundamental y todos los armónicos.

Por el momento puedo resolver el punto 1 de manera muy cruda al encontrar el pico más grande en el periodograma, luego multiplicarlo por donde es el número estimado de picos armónicos, pero estoy seguro de que debe haber Ser un mejor método que este truco.$1, 2, ..., N$$N$

Con respecto al punto 2, me encontré con esta pregunta que menciona el uso de un filtro de peine para eliminar un fundamental y todos sus armónicos. ¿Es este el mejor método para usar? Una consideración importante es que tendré que aplicar el filtro a cada píxel de series de tiempo en una gran matriz, por lo que un método computacionalmente eficiente sería altamente deseable.

Datos de ejemplo

• He subido una serie de tiempo 1D de valores de píxeles medios de ejemplo en .csvformato aquí .
• También agregué una muestra de 1000 cuadros (diezmado espacialmente 2: 1) en un .matarchivo comprimido aquí

Su método no es tan malo para un primer intento.

Sin embargo, el siguiente método tiende a funcionar mejor:

1. Buscar máximos locales
2. Verifique los máximos que están cerca (espacio de 2 o 3 compartimientos) y combínelos
3. Crea algunas hipótesis sobre la frecuencia fundamental. Actualmente supone que el pico más alto es la frecuencia fundamental, que es una hipótesis. También debe verificar la posibilidad de que el pico más alto sea el primer armónico, es decir, que haya un pico más pequeño a la mitad de la frecuencia del pico principal. También puede tener otros casos a considerar, utilizando su conocimiento del problema en cuestión (¿Interferencia? ¿Latidos irregulares?).
4. Suponiendo cada una de esas hipótesis, encuentre la frecuencia del suelo ajustando la parábola a cada pico armónico. Cada pico producirá una estimación ligeramente diferente debido al ruido, pero estos errores no están correlacionados y son promedio. Una de las hipótesis conducirá a ajustes mucho mejores, elija la frecuencia de suelo predicha por esta hipótesis.
5. Usando la frecuencia del suelo que encontró en el paso 4 como un hecho, vuelva a colocar la parábola alrededor de cada pico para estimar la altura del pico. Tenga en cuenta que los picos probablemente caen entre dos contenedores.
6. Ahora tiene la ubicación y la fuerza de lo fundamental y sus armónicos, pero no la fase. Probablemente sea más fácil encontrar la fase del fundamental, restar eso, encontrar la fase del primer armónico, etc.

La razón principal por la que esto funciona mejor es el paso 4. Cualquier hipótesis incorrecta para una frecuencia de tierra fallará horriblemente cuando intente ajustar los picos alrededor de la ubicación de los armónicos predichos. Digamos que tienes un pico a 2 Hz. Esta podría ser la frecuencia del suelo o un primer armónico. Cuando prueba la "primera hipótesis armónica", es decir, si la frecuencia de tierra es realmente 1 Hz, ajusta la parábola a los datos alrededor de 1,2,3,4,5 ... Hz. Si esa hipótesis es incorrecta, obtienes basura alrededor de 1,3,5 Hz. Si es correcto, puede encontrar picos cerca de 1.1 Hz, 2.2, 3.3, 4.4 y 5.6, lo que sugiere que la frecuencia de tierra real es 1.12 Hz.

Estás buscando sustracción espectral iterativa . Aquí hay información general de un libro de Análisis de contenido de Alexander Lerch.

Recomendaría una autocorrelación en su periodograma. Puede construir filtros de muesca o pico negativo basados ​​en múltiplos de la frecuencia fundamental que produce la correlación.

Este código me ayudó a crear una buena gráfica de autocorrelación (cepstrum es otro buen método cuando los armónicos tienen más poder que el fundamental) http://note.sonots.com/SciSoftware/Pitch.html

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

Puede tomar un promedio local en el DFT para encontrar los latidos del corazón. Si cierto punto en este grupo es mayor thresholdy el mayor valor en algún rango a su alrededor, entonces es un latido o un armónico.

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

Podrías usar los índices de latidos ubicados previamente y unirlos linealmente. Si eso no es demasiado crudo para ti.

¿Cuál es el propósito de esta parte de la operación? ¿Es este un análisis médico que requiere la ausencia del latido del corazón para buscar otros (tal vez) datos ocultos? ¿O es esto solo por estética?

Me imagino que un filtro de peine será algo incontrolable. Se requerirá un ajuste cuidadoso del parámetro de retroalimentación. Podría ser viable.

EDITAR: ¿Requiere un algoritmo de detección de tono? Escribí uno de estos hace unos años ignorando los términos (cercanos a) DC y encontrando el primer máximo local. Luego, usando los contenedores a ambos lados, podría usar una interpolación cuadrática y encontrar la ubicación de la frecuencia máxima con mayor precisión que la resolución permitida por las frecuencias de los contenedores.

¿Sería más preciso utilizar armónicos? Probablemente solo agregaría una oportunidad para errores. Aunque puede intentar obtener un valor para lo fundamental con este método y encontrar el máximo en doble, triple, etc. y luego usar un método similar al anterior:

1. encuentre el máximo local 2. interpolado cuadrático para encontrar el valor inter-bin para la frecuencia.

Si vetas el paso cuadrático y solo tomas el máximo local, no obtendrás una frecuencia precisa y subir los armónicos ayudaría.