Filtrar una señal digital en línea en tiempo real usando python

7

Actualmente estoy tratando de aplicar un filtro de paso de banda a una señal en tiempo real. Hay muestras que llegan con una frecuencia de muestreo constante y me gustaría calcular la señal filtrada de paso de banda correspondiente.

Cuál sería la mejor forma de hacer esto? ¿Tengo que filtrar la totalidad (o al menos un poco) de la señal cada vez que ingresan algunas muestras nuevas o hay alguna forma (como el DFT deslizante) en la que es posible determinar de manera eficiente la nueva parte del filtrado? ¿señal?

Me gustaría usar un filtro de butterworth (para el análisis fuera de línea, actualmente estoy usando la mantequilla y el filtro de scipy). Sé que esta función puede devolver un retraso de filtro, pero no sé cómo usarla para obtener una señal constante.

BStadlbauer
fuente

Respuestas:

1

La mecánica fundamental de realizar el procesamiento de auido digital en tiempo real bajo una plataforma de PC de propósito general basada en X-windows se basa en el uso de una familia de arquitecturas de doble búfer.

En esta arquitectura, el sonido que llega a través de una entrada de micrófono / línea se convierte primero en muestras a través de la tarjeta de sonido ADC y luego se llena en un búfer de entrada a la frecuencia de muestreo seleccionada por el usuario Fs. Cuando este búfer está lleno, primero el hardware de la tarjeta de sonido notifica al sistema operativo y luego el sistema operativo notifica a su programa. Y su programa puede acceder al bloque y puede comenzar el procesamiento de muestras en ese bloque.

Sin embargo, al mismo tiempo que está ocupado con el bloque actual, su programa ya ha proporcionado otro (el segundo) búfer para que la tarjeta de audio llene esas muestras que están llegando mientras procesa el búfer previamente llenado. Cuando este búfer disponible actualmente se procesa por completo, debe comenzar a procesar el siguiente búfer inmediatamente sin demoras, lo cual es una necesidad fundamental para la reproducción de audio sin clics. De esta forma de doble búfer, tiene la posibilidad de crear un procesamiento de audio fluido sin fallas ni grietas.

Además, ya sea que realice un filtrado basado en FIR o IIR, puede filtrar todo el búfer de una vez como en el caso de FIR, o ir recursivamente muestra por muestra para un caso de IIR.

El tamaño del búfer es importante para el retraso inicial del procesamiento. Por lo tanto, si lo toma demasiado, debe esperar hasta que se llenen los dos búferes antes de emitir cualquier cosa. Por otro lado, si toma las memorias intermedias demasiado cortas, el sistema se verá abrumado por las interrupciones entrantes.

Una opción óptima es entre 128 y 1024 muestras. Estas longitudes de memoria intermedia son apropiadas para el procesamiento posterior de tipo FFT. También se puede aumentar el número de memorias intermedias para obtener un rendimiento más sólido en condiciones de carga del sistema variables. Pero se requieren al menos dos amortiguadores.

Grasa32
fuente
2
Aunque estoy procesando señales EEG, puedo aplicar esto perfectamente, ¡gracias!
BStadlbauer
1
Por cierto, esa es casi exactamente la descripción de la arquitectura de búfer en cascada de GNU Radio.
Marcus Müller
1
Me dirigí a su publicación en la extensión de mi respuesta, @ Fat32; Espero que les guste :)
Marcus Müller
@ MarcusMüller; Gracias por la cooperación. Lo aprecio;)
Fat32
5

¿Tengo que filtrar la totalidad (o al menos un poco) de la señal cada vez que ingresan algunas muestras nuevas o hay alguna forma (como el DFT deslizante) en la que es posible determinar de manera eficiente la nueva parte del filtrado? ¿señal?

Los filtros digitales no funcionan así, básicamente, los FIR o IIR clásicos pueden funcionar en cada muestra nueva . Realmente deberías leer qué son estos filtros y cómo las personas los modelan.

Me gustaría usar un filtro butterworth

Bueno, hay muchas implementaciones de eso por ahí,

Actualmente estoy usando scipy's butter and lfilter

de los cuales ya conoces uno!

Ahora, un filtro de Butterworth es algo recursivo, por lo que para calcular la siguiente parte de la señal muestreada, necesitará el último estado. Ese es exactamente el "estado de retraso del filtro zi" que lfilterregresa y puede tomar la siguiente llamada como ziparámetro.

pero no sé cómo usarlo para obtener una señal constante.

Creo que te refieres a "lograr un filtrado continuo".

Ahora, dicho esto, el punto es que te estás preparando para escribir tu propia arquitectura de transmisión. Yo no haría eso. Use un marco existente. Por ejemplo, está GNU Radio, que le permite definir gráficos de flujo de procesamiento de señales en Python, y también es inherentemente multiproceso, utiliza implementaciones de algoritmos altamente optimizados, tiene muchas facilidades de entrada y salida, y viene con una gran biblioteca de bloques de procesamiento de señales , que se puede escribir en Python o C ++, si es necesario hacerlo.

Por ejemplo, un diagrama de flujo que toma muestras de una tarjeta de sonido, las filtra y las escribe en un archivo es:

#!/usr/bin/env python2
# -*- coding: utf-8 -*-
##################################################
# GNU Radio Python Flow Graph
# Title: Butterworth Test
# Generated: Mon Feb  8 16:17:18 2016
##################################################

from gnuradio import audio
from gnuradio import blocks
from gnuradio import eng_notation
from gnuradio import filter
from gnuradio import gr
from gnuradio.eng_option import eng_option
from gnuradio.filter import firdes
from optparse import OptionParser


class butterworth_test(gr.top_block):

    def __init__(self):
        gr.top_block.__init__(self, "Butterworth Test")

        ##################################################
        # Variables
        ##################################################
        self.samp_rate = samp_rate = 48000

        ##################################################
        # Blocks
        ##################################################
        # taps from scipy.butter!
        self.iir_filter_xxx_0 = filter.iir_filter_ffd(([1.0952627450621233e-05, 0.00013143152940745496, 0.0007228734117410033, 0.0024095780391366808, 0.005421550588057537, 0.008674480940892064, 0.010120227764374086, 0.008674480940892081, 0.005421550588057554, 0.0024095780391366955, 0.0007228734117410089, 0.00013143152940745594, 1.0952627450621367e-05]), ([1.0, -4.4363862740719835, 10.215121830052535, -15.374408118154847, 16.57333784740102, -13.325056987818655, 8.133543488903097, -3.77641064765334, 1.3181452681671835, -0.3361758629961047, 0.05930166356243964, -0.0064815521348275, 0.00033130678123743994]), False)
        self.blocks_file_sink_0 = blocks.file_sink(gr.sizeof_float*1, "", False)
        self.blocks_file_sink_0.set_unbuffered(False)
        self.audio_source_0 = audio.source(samp_rate, "", True)

        ##################################################
        # Connections
        ##################################################
        self.connect((self.audio_source_0, 0), (self.iir_filter_xxx_0, 0))    
        self.connect((self.iir_filter_xxx_0, 0), (self.blocks_file_sink_0, 0))    

def main(top_block_cls=butterworth_test, options=None):

    tb = top_block_cls()
    tb.start()
    try:
        raw_input('Press Enter to quit: ')
    except EOFError:
        pass
    tb.stop()
    tb.wait()


if __name__ == '__main__':
    main()

Tenga en cuenta que este código se generó automáticamente a partir de un gráfico de flujo gráfico en el que acabo de hacer clic usando el gnuradio-companionprograma:

diagrama de flujo diseñado en GRC

Si desea saber más sobre cómo implementar gráficos de flujo de procesamiento de señales en Python, vaya a los Tutoriales guiados por radio de GNU .

EDITAR : ¡Me gustó mucho la respuesta de @ Fat32! Lo que él describe como una arquitectura de doble buffer es bastante similar a lo que hace GNU Radio:

Un bloque ascendente produce muestras en fragmentos de muestra de tamaños arbitrarios, las escribe en el búfer de anillo de salida (que se representa como una flecha en la imagen de arriba) y notifica a sus bloques descendentes que hay nuevos datos.

El bloque aguas abajo recibe una notificación, comprueba si hay suficiente espacio en su búfer de salida para procesar las muestras que están en su búfer de anillo de entrada (que es lo mismo que el búfer de salida del bloque aguas arriba), los procesa. Cuando finaliza, informa a los bloques ascendentes que ha utilizado el búfer de anillo de entrada (que los bloques ascendentes pueden reutilizar como salida), y los bloques descendentes sobre las nuevas muestras disponibles.

Ahora, con GNU Radio siendo multiproceso, el bloque ascendente podría estar produciendo muestras nuevamente; En una aplicación normal de GNU Radio, casi todos los bloques están "activos" simultáneamente y las cosas escalan bastante bien en máquinas con múltiples CPU.

Por lo tanto, el trabajo principal de GNU Radio es brindarle esta infraestructura de búfer, la notificación y el mantenimiento de hilos, el bloque de procesamiento de señal claro API y algo para definir cómo está conectado todo, para que no tenga que escribir lo que Fat32 describe en su / su publícate! Tenga en cuenta que hacer una clasificación de flujo de muestra no es tan fácil de hacer correctamente, y GNU Radio le quita la dureza y le permite concentrarse en lo que quiere hacer: DSP.

Marcus Müller
fuente
¡Gracias! He examinado GNU Radio, pero como procesaré una señal de EEG, tendría que construir mi propio módulo para usar el gráfico de flujo porque cada muestra tiene una marca de tiempo que debe ser rastreable durante todo el proceso de filtrado.
BStadlbauer
No necesitas un bloqueo para eso. Las muestras son consecutivas, y con una frecuencia de muestreo fija, el tiempo está directamente disponible como la frecuencia de muestreo actual
Marcus Müller
@ MarcusMüller; Esta arquitectura de GNU Radio, que usted describe, es realmente lo que la filosofía moderna ofrece para su uso. Flexibilidad, facilidad de codificación y, lo más importante, poder concentrarse en cuál es su objetivo principal (procesamiento DSP), en lugar de cómo lograrlo utilizando detalles intrincados de bajo nivel (como lo que sucede cuando intenta implementar la técnica de doble buffer usando llamadas API Win32 !)
Fat32
@ Fat32 deberíamos atenuar el estilo de marketing pero: Sí, de hecho, y todo esto junto con ofrecer una arquitectura de búfer de anillo de copia cero eficiente y una extensa biblioteca de bloques que usan código optimizado a mano en x86, es MMX, SSE, SSE2, extensiones AVX y NEON de ARM donde corresponda usando la Biblioteca de Kernels Vector Optimizada VOLK :)
Marcus Müller
@ MarcusMüller Ayer tuve algo de tiempo libre, así que eché un vistazo más profundo a GNU Radio y me parecería bastante útil, solo que no entiendo cómo puedo "manualmente" alimentar muestras a los bloques sp, porque (todos) los turorials son hecho para fuentes de audio provenientes de un dispositivo de hardware. ¿No sabes dónde puedo encontrar un turial, etc. para empujar muestras manualmente a la cadena? PS, las muestras de EEG provienen de una capa especial (LSL - Lab Streaming Layer) y tienen alrededor de 64 canales por muestra
BStadlbauer