¿Como ha ido el verano? Nosotros trabajando bastante, aunque suene absurdo. Este mes de Agosto hemos actualizado poco porque hemos estado bastante liados, y entre todo el trabajo una de las cosas que hemos tenido que hacer es implementar filtros FIR en una FPGA, así que aprovechando que lo tenemos fresco (verano, fresco, jeje), hemos preparado este tutorial paso a paso sobre como diseñar, crear en VHDL y simular un filtro.

Antes de nada, para seguir este tutorial, aparte de las herramientas de siempre, Modelsim SE y Vivado, necesitas Python. ¿Por qué hace falta Python? Pues porque lo vamos a utilizar para crear las muestras que utilizaremos para simular nuestro filtro en VHDL y también para visualizar los resultados. Yo uso Python porque es gratis, pero podéis usar Matlab, Octave o lo que os de la gana. En concreto para instalar Python utilizo Anaconda, que te permite crear un Notebook desde donde ejecutar Python a través de Jupiter. No os preocupéis que es super fácil, te instalas Anaconda y cuando se abra seleccionas JupyterLab entre las cosas que salen y ya estas listo para empezar.

Pero antes de empezar a escribir código, habrá mucho lo prometo, vamos a definir que es un filtro FIR para los despistados.

¿Que es un filtro FIR?

Un filtro FIR (Finite Impulse Response) es una estructura de filtrado utilizada en procesado digital de señal, que se caracteriza por tener una respuesta al impulso finita en el tiempo, es decir que si dejas de meterle entradas, en algún momento la salida se convierte en cero. En un filtro IIR (Infinite impulse Response) la salida la tienes que cortar porque no se acaba nunca. Los filtros FIR se utilizan mucho por sus propiedades que no voy a ponerme a explicar (fase lineal, mucho mas fáciles de implementar, ..). Esto es un blog de VHDL y FPGA no de matemáticas y si estáis leyendo como se hace un filtro en VHDL, imagino que es porque o bien os lo han mandado de deberes en la universidad o tenéis que hacer uno en el trabajo y ya sabéis lo que es.

Para lo que a nosotros nos importa, lo que necesitamos es saber que el filtro puede ser paso bajo, paso alto o paso banda, y que tiene unas frecuencias de corte que es donde esta el limite de las frecuencias que deja pasar o eliminar.

Diseñando el filtro a implementar

Lo primero que necesitamos es un filtro. Para trabajar vamos a crear un filtro sencillo paso bajo. La señal que vamos a filtrar tiene una frecuencia máxima de 300Hz. Para redondear haremos un filtro entre 0 y 1000Hz con una frecuencia de corte de 100 Hz. Para ello vamos a ir a la página fiiir.com.

Los parametros del filtro serán:

  • Sampling freq 2000. Es el doble que nuestra frecuencia máxima de 1000Hz, según el teorema de Nyquist.
  • Cutoff frequency 100Hz. El filtro solo dejara pasar frecuencias hasta 100Hz y mayores las eliminaran.
  • Transition Bandwidth. Hemos escogido 300Hz. Es la velocidad con la que la pendiente del filtro cae cuando llegamos a la frecuencia de corte. Si la ponemos mas pequeña el filtro será mas restrictivo, pero el número de coeficientes crecerá bastante, y con más coeficientes necesitaremos mas hardware en la FPGA. Hay que elegirla con cuidado según las características de la señal de entrada.
  • Windows type. La restangular que hemos escogido significa no aplicar nada. Tenéis otras opciones que hacen que el filtro sera mucho mejor, pero de nuevo incrementa mucho el número de coeficientes y queremos un filtro sencillo. Jugad con ello para verlo.

Le dais a calcular y ya tenéis los 7 coeficientes del filtro que vamos a usar.

Diseño del filtro a implementar

Simulación en Python

Una vez tenemos los coeficientes vamos a probar el filtro y ver si funciona. Para ello abrimos el notebook de Python y creamos la señal de entrada y la pintamos.

import numpy as np
import matplotlib.pyplot as plt

#Frecuencia de muestreo 2000Hz
#Filtro respuesta de 0 a 1000
#Corte del filtro a 100Hz

#Tomados de fiiir.com
taps = np.array([
    0.130951924144877080,
    0.142713213023077851,
    0.150057552393435101,
    0.152554620877219965,
    0.150057552393435101,
    0.142713213023077851,
    0.130951924144877080,
])

t = np.linspace(0,1.0,2001)

sin_50Hz = np.sin(2*np.pi*50*t)
sin_300Hz = np.sin(2*np.pi*300*t)

original_signal = sin_50Hz+sin_300Hz
#Reescale to 0-1
original_signal = original_signal /2

plt.plot(t[:500],original_signal[:500])

Este código nos genera una señal que contiene una señal sinusoidal de 50Hz y otra de 300Hz sumadas. Cuando apliquemos el filtro con corte en 100Hz lo que haremos será cargarnos la señal de 300Hz.

Señal original

Vamos a filtrar la señal y ver que sale. Para ello tomaremos la entrada y haremos la convolución con el filtro. Como se ve el resultado es que hemos eliminado la señal de 300Hz y los quedamos solo con la de 50Hz. Parece que el filtro funciona.

filtered_signal = np.convolve(original_signal,taps)

plt.plot(t[:500],filtered_signal[:500])
Señal filtrada

Ahora ya podríamos pasar a VHDL a diseñar el filtro, pero tenemos un problema. Estos coeficientes del filtro y entradas son números en coma flotante y aunque es posible hacer el filtro en coma flotante en VHDL, no es muy buena idea ya que ocupará muchísimo más y será mucho mas lento. Así que ahora tenemos que tomar los coeficientes y la señal de entrada y pasarlos a punto fijo.

Para pasar a punto fijo, lo primero es decidir cuantos bits vamos a utilizar. En ese caso tenemos valores entre -1 y 1, así que podemos escoger por ejemplo 16 bits para representar nuestros valores en punto fijo.

El siguiente código Python, toma nuestros coeficientes y señal de entrada y los convierte a punto fijo multiplicando por 2^15 y luego haciendo la conversión a un valor int16 (entero con signo).

Después lo que tenemos es ya directamente el bucle que realiza el filtrado. En el vamos recorriendo todas las muestras de la señal de entrada y aplicamos el filtro que lo que hace es multiplicar las muestras anteriores de la entrada por cada uno de los coeficientes y sumarlos. Cosas a tener en cuenta:

  • Hay que saturar el resultado, para eso son los dos if. Si el valor se sale de rango lo cortamos para que no cambie de signo y distorsione la salida completamente.
  • Al multiplicar y acumular, si la entrada es de 16 bits necesitamos una salida de 32 bits Si os fijáis en en bucle hacemos la extensión. Esto habrá que hacerlo también en VHDL.
  • La salida la volvemos a convertir a float dividiento por el factor Q, en este caso 30, ya que la salida es un valor en punto fijo con 2 bits de parte entera y 30 de fraccionaria.

Si dubujamos el resultado del filtrado obtendremos exactamente lo mismo que en el caso de la implementación con flotantes y la convolución. Nuestro algoritmo en punto fijo funciona y ya podemos pasarlo a VHDL.

Q = 2**15
taps_Q = taps * Q
original_signal_Q = original_signal * Q

taps_fixed = np.int16(taps_Q)
#print('\n'.join([hex(i) for i in taps_fixed]))

original_signal_fixed = np.int16(original_signal_Q)
#for i in range(500):
#    print(original_signal_fixed[i])

#plt.plot(t[:100],original_signal_fixed[:100])

f_fixed = np.zeros(len(original_signal_fixed)+len(taps_fixed))

for j in range(len(original_signal_fixed+len(taps_fixed))):
    tmp = 0;
    for i in range(len(taps_fixed)):
        if((j-i)>0):        
            tmp = np.int32(tmp) + np.int32(taps_fixed[i])*np.int32(original_signal_fixed[j-i])
        if(tmp>0x3fffffff):
            tmp=0x3fffffff
        if(tmp<-0x40000000):
            tmp=-0x40000000      
    f_fixed[j] = tmp


f_float=f_fixed/(2**30)    

plt.plot(t[:500],f_float[:500])

Implementación en VHDL

Ya tenemos nuestro algoritmo y los coeficientes y señal de entrada para simularlo. Para generar las entradas solo tenéis que descomentar los print del código anterior en Python. Con esos valores vamos a trabajar.

Hay montones de formas de implementar el filtro FIR en VHDL dependiendo del rendimiento que necesitemos y el área que queramos utilizar. Nosotros vamos a hacer el diseño más sencillo para que podáis entenderlo y trabajar con el. Este diseño es totalmente segmentado, es decir, se puede introducir una muestra en cada ciclo de reloj y nos dará una salida también cada ciclo de reloj. Al tener 7 coeficientes, cada entrada tendrá que atravesar 7 registros hasta la salida. Es una implementación grande pero muy potente.

library ieee;
use ieee.std_logic_1164.all;
use ieee.numeric_std.all;

entity fir_filter is
port(
  clk_i : in std_logic;
  reset_i : in std_logic;
  x       : in signed(15 downto 0);
  y       : out signed(31 downto 0)
);
end entity;


architecture behav of fir_filter is

  type fixed_t is array (0 to 6) of signed(15 downto 0);
  signal taps : fixed_t := (X"10c3",X"1244",X"1335",X"1386",X"1335",X"1244",X"10c3");

  signal delay : fixed_t;

begin

delay_proc : process(clk_i, reset_i)
begin
  if(reset_i='1') then
    for i in 0 to 6 loop
      delay(i) <= (others=>'0');
    end loop;
  elsif(rising_edge(clk_i)) then
    delay(0) <= x;
    for i in 1 to 6 loop
      delay(i) <= delay(i-1);
    end loop;
  end if;
end process;

calc_proc : process(clk_i, reset_i)
  variable acum : signed(31 downto 0);
begin
  if(reset_i='1') then
    y <= (others=>'0');
    acum := (others=>'0');
  elsif (rising_edge(clk_i)) then
    acum:=(others=>'0');
    for i in 0 to 6 loop
      acum := acum + taps(i)*delay(i);
    end loop;
    y <= acum;
  end if;
end process;

end architecture;


library ieee;
use ieee.std_logic_1164.all;
use ieee.numeric_std.all;
use std.textio.all;

entity fir_tb is
end entity;


architecture behav of fir_tb is

  signal clk_i   : std_logic :='0';
  signal reset_i : std_logic;
  signal x       : signed(15 downto 0);
  signal y       : signed(31 downto 0);

  file samples : TEXT open READ_MODE is "samples.txt";
  file filtered : TEXT open WRITE_MODE is "output.txt";
begin

UUT : entity work.fir_filter(behav) port map(clk_i,reset_i,x,y);

clk_i   <= not(clk_i) after 10 ns;
reset_i <= '0', '1' after 10 ns, '0' after 20 ns;

process(clk_i)
  variable sample_line : LINE;
  variable x_int : integer;

begin
  if(rising_edge(clk_i)) then
    readline(samples,sample_line);
    read(sample_line,x_int);
    x <= to_signed(x_int,16);
  end if; 
end process;


process(clk_i)
  variable output_line : LINE;
  variable x_int : integer;

begin
  if(falling_edge(clk_i)) then
    write(output_line,to_integer(y));
    writeline(filtered,output_line);

  end if; 
end process;

end architecture;

Como veis cargamos los 7 coeficientes que hemos generado en Python en punto fijo en un array de registros. Luego hay dos procesos, uno que registra las entradas y las almacena retrasadas y otro proceso que multiplica esas entradas retrasadas con cada uno de los coeficientes. Es muy sencillo.

El testbench que adjuntamos toma las entradas de un fichero que se llama samples.txt donde hemos copiado los valores que hemos obtenido de imprimir la señal sinusoidal en punto fijo generada en Python. Leemos cada ciclo un valor y lo metemos en la entrada del filtro. La salida se guarda en un fichero de salida para poder cargarlo en Python y dibujarlo.

Si ejecutamos en Modelsim el sistema las primeras muestras del fichero de salida son:

0
0
0
67870747
162547987
227950000
239019507
221835098
231066891 

Que coinciden, con las que obtenemos del print de la señal de salida en punto fijo que hacemos en Python, con un pequeño error de redondeo, ya que en nuestro VHDL no redondeamos.

 [0.00000000e+00 6.78707470e+07 1.62547987e+08 2.27950000e+08  2.39019507e+08 2.21835098e+08 2.31066891e+08 3.03593593e+08  3.66540735e+08 4.16872945e+08] 

Como veis el filtro funciona bien. Asi que vamos a irnos a Vivado e implementarlo, a ver que sale. El resultado de la síntesis es el siguiente:

Como se puede ver el resultado es grande. Tenemos montones de registros, lógico si pensamos que tenemos 7 registros de 16 bits donde almacenamos los datos de entrada retrasados. Además podemos ver esos 7 cuadrados que son los multiplicadores/acumuladores que usamos para multiplicar las entradas con sus coeficientes. El resultado es el que esperabamos.

Esta es la implementación directa del filtro. Existen muchas otras formas de hacerlo y arquitecturas más óptimas. Si queréis hablar sobre ellas los comentarios están abiertos.