diff --git a/invisible_cities/cities/components.py b/invisible_cities/cities/components.py
index 3d18490bca..8a55f4c8fb 100644
--- a/invisible_cities/cities/components.py
+++ b/invisible_cities/cities/components.py
@@ -125,7 +125,7 @@ def create_timestamp(rate: float) -> float:
 
     Returns
     -------
-    Function to calculate timestamp for the given rate with 
+    Function to calculate timestamp for the given rate with
     event_number as parameter.
     """
 
@@ -142,12 +142,12 @@ def create_timestamp(rate: float) -> float:
     def create_timestamp_(event_number: Union[int, float]) -> float:
         """
         Calculates timestamp for a given Event Number and Rate.
-    
+
         Parameters
         ----------
         event_number : Union[int, float]
                        ID value of the current event.
-    
+
         Returns
         -------
         Calculated timestamp : float
@@ -326,6 +326,33 @@ def deconv_pmt(RWF):
     return deconv_pmt
 
 
+def deconv_pmt_fpga(dbfile     : str
+                   ,run_number : int
+                   ,selection  : List[int] = None
+                   ) -> Callable:
+    '''
+    Apply deconvolution as done in the FPGA.
+
+    Parameters
+    ----------
+    dbfile     : Database to be used
+    run_number : Run number of the database
+    selection  : List of PMT IDs to apply deconvolution to.
+
+    Returns
+    ----------
+    deconv_pmt : Function that will apply the deconvolution.
+    '''
+    DataPMT    = load_db.DataPMT(dbfile, run_number = run_number)
+    pmt_active = np.nonzero(DataPMT.Active.values)[0].tolist() if selection is None else selection
+    coeff_c    = DataPMT.coeff_c  .values.astype(np.double)[pmt_active]
+    coeff_blr  = DataPMT.coeff_blr.values.astype(np.double)[pmt_active]
+    def deconv_pmt(RWF):
+        return zip(*map(blr.deconvolve_signal_fpga, RWF[pmt_active],
+                        coeff_c                   , coeff_blr      ))
+    return deconv_pmt
+
+
 def get_run_number(h5in):
     if   "runInfo" in h5in.root.Run: return h5in.root.Run.runInfo[0]['run_number']
     elif "RunInfo" in h5in.root.Run: return h5in.root.Run.RunInfo[0]['run_number']
diff --git a/invisible_cities/cities/valdrada.py b/invisible_cities/cities/valdrada.py
new file mode 100644
index 0000000000..10e8ff842f
--- /dev/null
+++ b/invisible_cities/cities/valdrada.py
@@ -0,0 +1,172 @@
+"""
+-----------------------------------------------------------------------
+                                 valdrada
+-----------------------------------------------------------------------
+
+This city finds simulates the trigger procedure at the FPGA over PMT waveforms.
+This includes a number of tasks:
+    - Use a moving average to obtain the baseline
+    - Truncate to integers the deconvolution output
+    - Find trigger candidates and compute their characteristics, taking into
+      account the unique circumstances derived from working on a continuous scheme
+      (for example, delay in signals)
+    - Evaluate coincidences between trigger candidates.
+
+A number of general configuration parameters are needed (input as a dict, trigger_config):
+- coincidence_window : Time window (in time bins) in which valid trigger coincidences are counted
+- discard_width      : Any trigger with width less than this parameter is discarded.
+- multipeak          : Dictionary with multipeak protection parameters, otherwise None.
+     - q_min      : Integrated charge threshold of a post-trigger peak to discard
+                    a trigger due to multipeak protection. In ADC counts.
+     - time_min   : Minimum width of a post-trigger peak to discard a trigger due to
+                    multipeak protection. In time bins.
+     - time_after : For how long is the multipeak protection evaluated after a
+                    valid trigger. In mus.
+
+A individual trigger configuration can be given per channel, through a dict
+with keys equal to PMT IDs, which marks the validity range of the peak
+characteristics:
+- q_min, q_max          : Range for the integrated charge of the peak (q_min < q < q_max).
+                          In ADC counts.
+- time_min, time_max    : Range for the peak width (time_min <= width < time_max).
+                          In time mus.
+- baseline_dev, amp_max : Range for peak height (baseline_dev < height < amp_max).
+                          In ADC counts.
+- pulse_valid_ext       : Time allowed for a pulse to go below baseline_dev. In ns.
+
+The result of the city is a dataframe containing the event ID and PMT ID of each
+trigger candidate. For each trigger candidate a number of parameters is computed:
+ - trigger_time : time bin at which the trigger candidate starts.
+ - q            : integrated ADC counts within the peak.
+ - width        : Length of the peak in time bins.
+ - height       : Maximum ADC values in a given time bien within the peak.
+ - baseline     : Value of the baseline at trigger_time.
+ - max_height   : Maximum (minimum due to PMT negative signals) height in the wvf.
+
+Additionally, a set of flags are assigned depending on wether the parameters are
+within range of the trigger configuration:
+ - valid_q    : If peak is within q range.
+ - valid_w    : If peak is within width range.
+ - valid_h    : If peak is within height range.
+ - valid_peak : Only if 'multipeak' is active, True if there isn't a post-trigger
+               candidate with the configuration parameters.
+ - valid_all  : boolean and of previous valid flags.
+
+ Finally, a series of coincidence-related values are also given:
+ - n_coinc       : Number of valid triggers within the coincidence_window,
+                   starting at the trigger trigger_time and including the trigger itself.
+                   -1 indicates no valid trigger (including the trigger itself).
+ - closest_ttime : Time difference to closest valid trigger.
+                   -1 if there are none aside from the trigger itself.
+ - closest_pmt   : PMT ID of the closest valid trigger.
+                   -1 if there are none aside from the trigger itself.
+"""
+
+import tables as tb
+import numpy  as np
+
+from .. reco                  import tbl_functions        as tbl
+from .. io  .run_and_event_io import run_and_event_writer
+from .. io  .trigger_io       import       trigger_writer, trigger_dst_writer
+
+from .. dataflow            import dataflow as fl
+from .. dataflow.dataflow   import push
+from .. dataflow.dataflow   import pipe
+from .. dataflow.dataflow   import sink
+
+from .  components import city
+from .  components import print_every
+from .  components import collect
+from .  components import copy_mc_info
+from .  components import deconv_pmt_fpga
+from .  components import WfType
+from .  components import wf_from_files
+from .  components import get_number_of_active_pmts
+
+from .. reco.trigger_functions import retrieve_trigger_information
+from .. reco.trigger_functions import get_trigger_candidates
+from .. reco.trigger_functions import check_trigger_coincidence
+
+from .. core       import system_of_units as units
+
+
+def check_empty_trigger(triggers) -> bool:
+    """
+    Filter for valdrada flow.
+    The flow stops if there are no candidate triggers/
+    """
+    return len(triggers) > 0
+
+@city
+def valdrada(files_in, file_out, compression, event_range, print_mod, detector_db, run_number,
+             trigger_config = dict(), channel_config = dict()):
+
+    # Change time units into number of waveform bins.
+    if trigger_config['multipeak'] is not None:
+        trigger_config['multipeak']['time_min'  ] /= 25*units.ns
+        trigger_config['multipeak']['time_after'] /= 25*units.ns
+    for k in channel_config.keys():
+        channel_config[k]['time_min']       /= 25*units.ns
+        channel_config[k]['time_max']       /= 25*units.ns
+        channel_config[k]['pulse_valid_ext'] = int(channel_config[k]['pulse_valid_ext']/25*units.ns)
+
+    #### Define data transformations
+    # Raw WaveForm to deconvolved waveforms
+    rwf_to_cwf           = fl.map(deconv_pmt_fpga(detector_db, run_number, list(channel_config.keys())),
+                                  args = "pmt",
+                                  out  = ("cwf", "baseline"))
+
+    # Extract all possible trigger candidates of each PMT
+    trigger_on_channels  = fl.map(retrieve_trigger_information(channel_config, trigger_config),
+                                  args = ("pmt", "cwf", "baseline", "event_number"),
+                                  out  = "triggers")
+
+    # Add coincidence between channels
+    check_coincidences   = fl.map(check_trigger_coincidence(trigger_config['coincidence_window']),
+                                  item = "triggers")
+
+    # Filter events with zero triggers
+    filter_empty_trigger = fl.map(check_empty_trigger,
+                                  args = "triggers",
+                                  out  = "empty_trigger")
+
+    event_count_in     = fl.spy_count()
+    event_count_out    = fl.spy_count()
+    events_no_triggers = fl.count_filter(bool, args = "empty_trigger")
+    evtnum_collect     = collect()
+
+    with tb.open_file(file_out, "w", filters = tbl.filters(compression)) as h5out:
+        # Define writers...
+        write_event_info_   = run_and_event_writer(h5out)
+        write_trigger_info_ = trigger_writer      (h5out, get_number_of_active_pmts(detector_db, run_number))
+        write_trigger_dst_  = trigger_dst_writer  (h5out)
+        # ... and make them sinks
+
+        write_event_info   = sink(write_event_info_  , args=( "run_number"  , "event_number"    , "timestamp"   ))
+        write_trigger_info = sink(write_trigger_info_, args=( "trigger_type", "trigger_channels"                ))
+        write_trigger_dst  = sink(write_trigger_dst_ , args=  "triggers"                                         )
+
+        result = push(source = wf_from_files(files_in, WfType.rwf),
+                      pipe   = pipe(fl.slice(*event_range, close_all=True),
+                                    print_every(print_mod),
+                                    event_count_in.spy,
+                                    rwf_to_cwf,
+                                    trigger_on_channels,
+                                    filter_empty_trigger,
+                                    events_no_triggers.filter,
+                                    check_coincidences,
+                                    event_count_out.spy,
+                                    fl.branch("event_number", evtnum_collect.sink),
+                                    fl.fork(write_trigger_dst ,
+                                            write_event_info  ,
+                                            write_trigger_info)),
+                      result = dict(events_in   = event_count_in    .future,
+                                    events_out  = event_count_out   .future,
+                                    evtnum_list = evtnum_collect    .future,
+                                    events_pass = events_no_triggers.future))
+
+        if run_number <= 0:
+            copy_mc_info(files_in, h5out, result.evtnum_list,
+                         detector_db, run_number)
+
+        return result
diff --git a/invisible_cities/cities/valdrada_test.py b/invisible_cities/cities/valdrada_test.py
new file mode 100644
index 0000000000..d6da771cbc
--- /dev/null
+++ b/invisible_cities/cities/valdrada_test.py
@@ -0,0 +1,53 @@
+import os
+import tables as tb
+
+from .  valdrada           import valdrada
+from .. core.testing_utils import assert_tables_equality
+
+
+def test_valdrada_contains_all_tables(trigger_config):
+    conf, PATH_OUT = trigger_config
+    valdrada(**conf)
+    with tb.open_file(PATH_OUT) as h5out:
+        assert "Run"             in h5out.root
+        assert "Run/events"      in h5out.root
+        assert "Run/runInfo"     in h5out.root
+        assert "Trigger"         in h5out.root
+        assert "Trigger/events"  in h5out.root
+        assert "Trigger/trigger" in h5out.root
+        assert "Trigger/DST"     in h5out.root
+
+
+def test_valdrada_exact_result_multipeak(ICDATADIR, trigger_config):
+    true_out          = os.path.join(ICDATADIR, "exact_result_multipeak_valdrada.h5")
+    conf, PATH_OUT   = trigger_config
+    valdrada(**conf)
+
+    tables = ("Trigger/events", "Trigger/trigger", "Trigger/DST",
+              "Run/events"      , "Run/runInfo")
+
+    with tb.open_file(true_out)  as true_output_file:
+        with tb.open_file(PATH_OUT) as      output_file:
+            for table in tables:
+                assert hasattr(output_file.root, table)
+                got      = getattr(     output_file.root, table)
+                expected = getattr(true_output_file.root, table)
+                assert_tables_equality(got, expected)
+
+
+def test_valdrada_exact_result(ICDATADIR, trigger_config):
+    true_out          = os.path.join(ICDATADIR, "exact_result_valdrada.h5")
+    conf, PATH_OUT    = trigger_config
+    conf['trigger_config']['multipeak'] = None
+    valdrada(**conf)
+
+    tables = ("Trigger/events", "Trigger/trigger", "Trigger/DST",
+              "Run/events"      , "Run/runInfo")
+
+    with tb.open_file(true_out)  as true_output_file:
+        with tb.open_file(PATH_OUT) as      output_file:
+            for table in tables:
+                assert hasattr(output_file.root, table)
+                got      = getattr(     output_file.root, table)
+                expected = getattr(true_output_file.root, table)
+                assert_tables_equality(got, expected)
diff --git a/invisible_cities/config/valdrada.conf b/invisible_cities/config/valdrada.conf
new file mode 100644
index 0000000000..cfd0424cf3
--- /dev/null
+++ b/invisible_cities/config/valdrada.conf
@@ -0,0 +1,19 @@
+files_in    = '$ICDIR/database/test_data/blr_fpga_examples.h5'
+file_out   = '$ICDIR/database/test_data/valdrada.h5'
+
+compression = 'ZLIB4'
+run_number  = 8093
+detector_db = 'new'
+event_range = all
+
+trigger_config = {'coincidence_window':64, 'discard_width':40,
+                  'multipeak':{'q_min':100000, 'time_min':2*mus, 'time_after':800*mus}}
+
+channel_config = {0: {'q_min'       :5000, 'q_max'   :50000,
+                      'time_min'    :2*mus, 'time_max': 40*mus,
+                      'baseline_dev':  10, 'amp_max' : 1000,
+                      'pulse_valid_ext':50*ns},
+                  2: {'q_min'       :5000, 'q_max'   :50000,
+                      'time_min'    :2*mus, 'time_max': 40*mus,
+                      'baseline_dev':  10, 'amp_max' : 1000,
+                      'pulse_valid_ext':50*ns}}
diff --git a/invisible_cities/conftest.py b/invisible_cities/conftest.py
index 32bfbd1aeb..fd81ed4040 100644
--- a/invisible_cities/conftest.py
+++ b/invisible_cities/conftest.py
@@ -50,6 +50,11 @@ def example_blr_wfs_filename(ICDATADIR):
     return os.path.join(ICDATADIR, "blr_examples.h5")
 
 
+@pytest.fixture(scope='session')
+def example_blr_fpga_wfs_filename(ICDATADIR):
+    return os.path.join(ICDATADIR, "blr_fpga_examples.h5")
+
+
 @pytest.fixture(scope  = 'session',
                 params = ['electrons_40keV_z250_MCRD.h5'])
 def electron_MCRD_file(request, ICDATADIR):
@@ -699,6 +704,41 @@ def deconvolution_config(ICDIR, ICDATADIR, PSFDIR, config_tmpdir):
 
     return conf, PATH_OUT
 
+
+@pytest.fixture(scope='function') # Needs to be function as the config dict is modified when running
+def trigger_config(ICDIR, ICDATADIR, config_tmpdir):
+    PATH_IN     = os.path.join(ICDATADIR    ,    "blr_fpga_examples.h5")
+    PATH_OUT    = os.path.join(config_tmpdir,     "trigger_city_out.h5")
+    nevt_req    = 10
+    conf        = dict(files_in       = PATH_IN ,
+                       file_out       = PATH_OUT,
+                       event_range    = nevt_req,
+                       compression    = 'ZLIB4',
+                       print_mod      = 1000,
+                       run_number     = 8093,
+                       trigger_config = dict(coincidence_window = 64
+                                            ,discard_width      = 40
+                                            ,multipeak          = dict(q_min      = 100000
+                                                                      ,time_min   = 2  *units.mus
+                                                                      ,time_after = 800*units.mus)),
+                       channel_config = {0 : dict(q_min           = 5000
+                                                 ,q_max           = 50000
+                                                 ,time_min        = 2 *units.mus
+                                                 ,time_max        = 40*units.mus
+                                                 ,baseline_dev    = 10
+                                                 ,amp_max         = 1000
+                                                 ,pulse_valid_ext = 50*units.ns)
+                                        ,2 : dict(q_min           = 5000
+                                                 ,q_max           = 50000
+                                                 ,time_min        = 2 *units.mus
+                                                 ,time_max        = 40*units.mus
+                                                 ,baseline_dev    = 10
+                                                 ,amp_max         = 1000
+                                                 ,pulse_valid_ext = 50*units.ns)})
+
+    return conf, PATH_OUT
+
+
 ## To make very slow tests only run with specific option
 def pytest_addoption(parser):
     parser.addoption(
diff --git a/invisible_cities/database/test_data/blr_fpga_examples.h5 b/invisible_cities/database/test_data/blr_fpga_examples.h5
new file mode 100755
index 0000000000..37be52aecb
--- /dev/null
+++ b/invisible_cities/database/test_data/blr_fpga_examples.h5
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:5acb7051d8dd3f7f0433baa616937089fd66d54d5ddd8bd470ad711d9be7cc0f
+size 15111797
diff --git a/invisible_cities/database/test_data/exact_result_multipeak_valdrada.h5 b/invisible_cities/database/test_data/exact_result_multipeak_valdrada.h5
new file mode 100644
index 0000000000..fe38b05357
--- /dev/null
+++ b/invisible_cities/database/test_data/exact_result_multipeak_valdrada.h5
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:0f11a28f2437e57f243594ee48aec241f821543f8d5fb3dd94545b4d6e2dbb5d
+size 21870
diff --git a/invisible_cities/database/test_data/exact_result_valdrada.h5 b/invisible_cities/database/test_data/exact_result_valdrada.h5
new file mode 100644
index 0000000000..54ce6b3c6c
--- /dev/null
+++ b/invisible_cities/database/test_data/exact_result_valdrada.h5
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:599e156701fc2e8b2cf35d3d421990feecb4448523afcae07c261eecadb769c5
+size 21863
diff --git a/invisible_cities/evm/nh5.py b/invisible_cities/evm/nh5.py
index 8d1be863fb..4c502f30d7 100644
--- a/invisible_cities/evm/nh5.py
+++ b/invisible_cities/evm/nh5.py
@@ -186,3 +186,22 @@ class VoxelsTable(tb.IsDescription):
 class EventPassedFilter(tb.IsDescription):
     event  = tb.Int32Col(pos=0)
     passed = tb. BoolCol(pos=1)
+
+
+class TriggerTable(tb.IsDescription):
+    event         = tb.UInt32Col(pos= 0)
+    pmt           = tb.UInt32Col(pos= 1)
+    trigger_time  = tb.UInt32Col(pos= 2)
+    q             = tb.UInt32Col(pos= 3)
+    width         = tb.UInt32Col(pos= 4)
+    height        = tb.UInt32Col(pos= 5)
+    valid_q       = tb.BoolCol  (pos= 6)
+    valid_w       = tb.BoolCol  (pos= 7)
+    valid_h       = tb.BoolCol  (pos= 8)
+    valid_peak    = tb.BoolCol  (pos= 9)
+    valid_all     = tb.BoolCol  (pos=10)
+    baseline      = tb.Int32Col (pos=11)
+    max_height    = tb.Int32Col (pos=12)
+    n_coinc       = tb.Int32Col (pos=13)
+    closest_ttime = tb.Int32Col (pos=14)
+    closest_pmt   = tb.Int32Col (pos=15)
diff --git a/invisible_cities/io/trigger_io.py b/invisible_cities/io/trigger_io.py
index 8f077e191c..236a2910b6 100644
--- a/invisible_cities/io/trigger_io.py
+++ b/invisible_cities/io/trigger_io.py
@@ -4,6 +4,7 @@
 
 from .. evm                import nh5     as table_formats
 from .. reco.tbl_functions import filters as tbl_filters
+from .. io  .table_io      import make_table
 
 
 def store_trigger(tables, trg_type, trg_channels):
@@ -41,3 +42,36 @@ def _make_tables(hdf5_file, n_sensors, compression="ZLIB4"):
     trg_tables = trg_type, trg_channels
 
     return trg_tables
+
+
+def trigger_dst_writer(hdf5_file, **kwargs):#{
+    trigger_table = make_table(hdf5_file,
+                               group       = "Trigger",
+                               name        = "DST"    ,
+                               fformat     = table_formats.TriggerTable,
+                               description = "Simulated trigger data",
+                               compression = 'ZLIB4')
+    def write_trigger(trigger_info):#{
+        row = trigger_table.row
+        row["event"        ] = trigger_info.event
+        row["pmt"          ] = trigger_info.pmt
+        row["trigger_time" ] = trigger_info.trigger_time
+        row["q"            ] = trigger_info.q
+        row["width"        ] = trigger_info.width
+        row["height"       ] = trigger_info.height
+        row["valid_q"      ] = trigger_info.valid_q
+        row["valid_w"      ] = trigger_info.valid_w
+        row["valid_h"      ] = trigger_info.valid_h
+        row["valid_peak"   ] = trigger_info.valid_peak
+        row["valid_all"    ] = trigger_info.valid_all
+        row["baseline"     ] = trigger_info.baseline
+        row["max_height"   ] = trigger_info.max_height
+        row["n_coinc"      ] = trigger_info.n_coinc
+        row["closest_ttime"] = trigger_info.closest_ttime
+        row["closest_pmt"  ] = trigger_info.closest_pmt
+        row.append()
+
+    def write_triggers(triggers):
+        for t in triggers: write_trigger(t)
+
+    return write_triggers
diff --git a/invisible_cities/io/trigger_io_test.py b/invisible_cities/io/trigger_io_test.py
new file mode 100644
index 0000000000..d30efdf11b
--- /dev/null
+++ b/invisible_cities/io/trigger_io_test.py
@@ -0,0 +1,82 @@
+import os
+
+import numpy  as np
+import tables as tb
+
+from numpy     .testing        import assert_allclose
+from hypothesis                import given
+from hypothesis.strategies     import integers
+from hypothesis.strategies     import lists
+from hypothesis.strategies     import booleans
+from hypothesis.strategies     import composite
+
+from ..core.testing_utils      import assert_dataframes_equal
+from ..io.dst_io               import load_dst
+from . trigger_io              import trigger_dst_writer
+from .. reco.trigger_functions import TriggerInfo
+
+
+@composite
+def trigger_value(draw):
+    uints = integers(min_value=0                     , max_value=np.iinfo(np.uint32).max)
+    ints  = integers(min_value=np.iinfo(np.int32).min, max_value=np.iinfo(np.int32).max)
+
+    # Event and PMT ID
+    evt_id = draw(uints)
+    pmt_id = draw(uints)
+
+    # Peak time, q, width and height
+    time   = draw(uints)
+    q      = draw(uints)
+    width  = draw(uints)
+    height = draw(uints)
+
+    # Trigger validity flags
+    valid_q    = draw(booleans())
+    valid_w    = draw(booleans())
+    valid_h    = draw(booleans())
+    valid_peak = draw(booleans())
+    valid_all  = all([valid_q, valid_w, valid_h, valid_peak])
+
+    # Baseline and wvf max. height
+    baseline   = draw(ints)
+    max_height = draw(ints)
+
+    # Number of coincidences, closest pmt and tbin
+    n_coinc      = draw(ints)
+    closest_time = draw(ints)
+    closest_pmt  = draw(ints)
+
+    trigger = [evt_id  , pmt_id    , time   , q           , width      , height
+              ,valid_q , valid_w   , valid_h, valid_peak  , valid_all
+              ,baseline, max_height, n_coinc, closest_time, closest_pmt        ]
+
+    return TriggerInfo(*trigger)
+
+
+@composite
+def trigger_values(draw):
+    size  = draw(integers(min_value=1, max_value=10))
+    return draw(lists(trigger_value(), min_size=size, max_size=size))
+
+
+@given(triggers=trigger_values())
+def test_trigger_dst_writer(config_tmpdir, triggers):
+    output_file = os.path.join(config_tmpdir, "test_trigger_dst.h5")
+
+    with tb.open_file(output_file, 'w') as h5out:
+        write = trigger_dst_writer(h5out)
+        write(triggers)
+
+    trigger_dst = load_dst(output_file, group='Trigger', node='DST')
+
+    col_names = ["event"     , "pmt"       , "trigger_time" , "q"
+                ,"width"     , "height"    , "valid_q"      , "valid_w"
+                ,"valid_h"   , "valid_peak", "valid_all"    , "baseline"
+                ,"max_height", "n_coinc"   , "closest_ttime", "closest_pmt"]
+
+    assert np.all([tname == cname for tname, cname in zip(trigger_dst.columns.values, col_names)])
+
+    # Check values stored
+    for i, row in trigger_dst.iterrows():
+        assert_allclose(triggers[i], list(row.values))
diff --git a/invisible_cities/reco/trigger_functions.py b/invisible_cities/reco/trigger_functions.py
new file mode 100644
index 0000000000..2e8c682010
--- /dev/null
+++ b/invisible_cities/reco/trigger_functions.py
@@ -0,0 +1,198 @@
+import numpy  as np
+
+from typing      import List
+from typing      import Dict
+from typing      import Callable
+
+from collections import namedtuple
+
+from .. core.core_functions import in_range
+
+TriggerInfo = namedtuple('TriggerInfo',
+                         'event pmt '
+                         'trigger_time q width height '
+                         'valid_q valid_w valid_h valid_peak valid_all '
+                         'baseline max_height n_coinc closest_ttime closest_pmt')
+
+def get_trigger_candidates(signal       : np.ndarray
+                          ,event_id     : int           , pmt_id         : int
+                          ,channel_conf : Dict[int, int], trigger_conf   : Dict
+                          ,baseline     : np.ndarray    , wvf_max_height : int
+                          ) -> List:
+    """
+    Calculates the trigger information of a given signal.
+
+    Parameters
+    ----------
+    signal         : Deconvolved waveform.
+    event_id       : Event ID associated to the waveform.
+    pmt_id         : PMT ID of the waveform's PMT.
+    channel_conf   : PMT specific trigger configuration.
+    trigger_conf   : General trigger configuration.
+    baseline       : Array with the baseline value at each time bin.
+    wvf_max_height : Maximum height of the waveform.
+
+    Returns
+    ----------
+    triggers       : List with all the trigger related information.
+    """
+    # Trigger parameters
+    qrange    = [channel_conf['q_min'          ], channel_conf['q_max'   ]]
+    wrange    = [channel_conf['time_min'       ], channel_conf['time_max']]
+    hrange    = [channel_conf['baseline_dev'   ], channel_conf['amp_max' ]]
+    pulse_ext =  channel_conf['pulse_valid_ext']
+
+    multipk       = trigger_conf['multipeak'    ] # Protection for later peaks
+    discard_width = trigger_conf['discard_width'] # Don't store events below certain width
+
+    # Multipeak protection parameters
+    if multipk:
+        multi_q =     multipk['q_min'     ]  # Minimum charge of later peaks to discard a trigger
+        multi_w = int(multipk['time_min'  ]) # Minimum time width of later peaks to discard a trigger
+        multi_t = int(multipk['time_after']) # Time window for a later peak to be considered
+
+    # Signal above threshold, allowing for pulse_ext counts below threshold
+    triggers   = []
+    indexes    = np.where(signal > hrange[0])[0]
+    if len(indexes) == 0: return triggers # Protection in case min. threshold is not surpassed
+    index_sel  = np.diff(indexes)>(pulse_ext+1)
+    candidates = np.split(indexes, np.where(index_sel)[0] + 1)
+
+    for p in candidates:
+        peak_slice = slice(p[0]-1, p[-1] + pulse_ext + 1) # DAQ gets pulse_ext counts after peaks; also, for Q it takes a bin before TODO Fix for N100 (remove -1)
+        peak       = signal[peak_slice]
+        if not (len(peak) > max(discard_width, 1)): continue
+        width      = len(peak) - 1           # TODO Remove -1 in N100
+        height     = peak[1:  ].max()        # TODO Full peak for N100
+        peak       = peak[ :-2][peak[:-2]>0] # Only counts above 0 are used for Q, last bin not considered TODO ONLY -1 for N100
+        charge     = peak.sum()
+        if len(peak)==0: continue
+        start_time = p[0]
+
+        # Trigger validity
+        trigger_w   = in_range(width , *wrange, left_closed=True , right_closed=False)
+        trigger_q   = in_range(charge, *qrange, left_closed=False, right_closed=False)
+        trigger_h   = height < hrange[1]
+        trigger_all = all([trigger_w, trigger_q, trigger_h])
+
+        # Baseline at the start of the pulse
+        baseline_t0 = baseline[start_time]
+
+        trigger     = TriggerInfo(event_id   , pmt_id
+                                 ,start_time , charge        , width    , height
+                                 ,trigger_q  , trigger_w     , trigger_h,   True, trigger_all  # By default multiflag is set as true
+                                 ,baseline_t0, wvf_max_height
+                                 ,-1         , -1            , -1) # Filled later in the coincidence module
+        triggers  .append(trigger)
+
+    # Extra-peaks protection
+    if multipk is not None:
+        # Evaluate if later peaks (must be over discard_width) are candidates for the multipeak protection
+        for i, t0 in zip(range(len(triggers[:-1])), triggers[:-1]):
+            new_val = t0.valid_peak
+            for t1 in triggers[i+1:]:
+                multi_candidate = ((t1.width        >=                              multi_w  ) &
+                                   (t1.q            >                               multi_q  ) &
+                                   (t1.trigger_time <= t0.trigger_time + t0.width + multi_t+2)) # Start of trigger + peak width + time window + FPGA delay (2)
+                if not multi_candidate:
+                    continue
+                # In case the later peak ends after the multipeak protection window ends
+                if (t1.trigger_time + t1.width) > (t0.trigger_time + t0.width + multi_t + 2): # Event partially contained in the multipeak window
+                    peak_slice = slice(t1.trigger_time - 1, t0.trigger_time + t0.width + multi_t + 2 + 1)
+                    peak       = signal[peak_slice]
+                    width      = len(peak) - 1 # TODO Remove -1 in N100
+                    charge     = peak[peak>0].sum()
+
+                    new_val    = (width >= multi_w) & (charge > multi_q)
+                else:
+                    new_val = False
+                if not new_val: break
+            triggers[i] = triggers[i]._replace(valid_peak=new_val, valid_all=all([triggers[i].valid_all, new_val]))
+
+    return triggers
+
+
+def retrieve_trigger_information(channel_config : Dict[int, int]
+                                ,trigger_config : Dict
+                                ) -> Callable:
+    """
+    Calculates the trigger information of deconvolved waveforms.
+
+    Parameters
+    ----------
+    channel_config : Channel-specific trigger configuration.
+    trigger_config : General trigger configuration.
+
+    Returns
+    ----------
+    trigger_on_channels : Function that returns each channel's trigger candidates.
+    """
+    def trigger_on_channels(rwfs         : np.ndarray
+                           ,cwfs         : np.ndarray
+                           ,baselines    : np.ndarray
+                           ,event_number : int):
+        triggers = []
+        for i, dict_items in enumerate(channel_config.items()):
+            pmt_id, pmt_conf = dict_items
+            baseline         = baselines[i]
+            signal           = cwfs[i] - baseline
+            max_height       = rwfs[pmt_id].min()
+            triggers.extend(get_trigger_candidates(signal
+                                                  ,event_number, pmt_id
+                                                  ,pmt_conf    , trigger_config
+                                                  ,baseline    , max_height    ))
+        return triggers
+    return trigger_on_channels
+
+
+def check_trigger_coincidence(coinc_window : int) -> Callable:
+    """
+    Checks if there is time coincidence between trigger candidates and modifies
+    the input list of trigger candidates accordingly.
+
+    Parameters
+    ----------
+    coinc_window   : Time window to considerate the coincidence condition as valid.
+
+    Returns
+    ----------
+    pmt_coincidence : Function that checks for coincidence and modififes triggers.
+    """
+
+    def pmt_coincidence(triggers : np.ndarray) -> np.ndarray:
+        valid_sel  = np.array([t.valid_all    for t in triggers])            # Use only valid triggers
+        times      = np.array([t.trigger_time for t in triggers])[valid_sel] # Times of valid triggers
+        ids        = np.array([t.pmt          for t in triggers])[valid_sel] # Ids   of valid triggers
+
+        # Order by time, index so we can use this to order both times and ids
+        order_idx  = np.lexsort((ids, times))
+        order_time = times[order_idx]
+        pmts       = np.full(ids.shape, -1)
+
+        # Calculate all time differences between each trigger and all later ones
+        time_diffs = [order_time[i+1:] - t for i, t in enumerate(order_time)]
+
+        # Since pmts are ordered by time, closest pmt will be the next one in the array.
+        pmts[:-1]  = ids[order_idx][1:]
+
+        # Shortest time is the first value of time_diff for each trigger
+        shortest   = np.array([t[0] if len(t)>0 else -1 for t in time_diffs])
+
+        # Count how many triggers are within the coincidince window of each trigger
+        n_coinc    = np.array([len(t[t<coinc_window])+1 for t in time_diffs])
+
+        # Undo the original ordering
+        old_order  = np.argsort(order_idx)
+
+        # Reorder the coincidence results
+        new_n_coinc      = n_coinc [old_order]
+        new_closest_time = shortest[old_order]
+        new_pmt_id       = pmts    [old_order]
+
+        # Update the triggers with the coincidence info
+        for j, i in enumerate(np.where(valid_sel)[0]):
+            triggers[i] = triggers[i]._replace(n_coinc      =new_n_coinc     [j],
+                                               closest_ttime=new_closest_time[j],
+                                               closest_pmt  =new_pmt_id      [j])
+        return triggers
+    return pmt_coincidence
diff --git a/invisible_cities/reco/trigger_functions_test.py b/invisible_cities/reco/trigger_functions_test.py
new file mode 100644
index 0000000000..35e68bd7f3
--- /dev/null
+++ b/invisible_cities/reco/trigger_functions_test.py
@@ -0,0 +1,331 @@
+import numpy as np
+
+from pytest                    import fixture
+from pytest                    import mark
+from hypothesis                import given
+from numpy     .testing        import assert_allclose
+from hypothesis.strategies     import integers
+from hypothesis.strategies     import lists
+from hypothesis.strategies     import booleans
+from hypothesis.strategies     import composite
+
+from .. reco.trigger_functions import get_trigger_candidates
+from .. reco.trigger_functions import retrieve_trigger_information
+from .. reco.trigger_functions import check_trigger_coincidence
+from .. reco.trigger_functions import TriggerInfo
+from .. core.core_functions    import in_range
+
+
+def channel_conf(**kwargs):
+    channel_dict = {'q_min'          : 5000, 'q_max'   :  50000,
+                    'time_min'       : 2000, 'time_max':  10000,
+                    'baseline_dev'   :    5, 'amp_max' :   1000,
+                    'pulse_valid_ext':   15}
+    channel_dict.update(kwargs)
+    return channel_dict
+
+
+def trigger_conf(**kwargs):
+    trigger_dict = {'coincidence_window':64, 'discard_width':40, 'multipeak':None}
+    trigger_dict.update(kwargs)
+    return trigger_dict
+
+
+@composite
+def double_square_wf(draw):
+    # Fake signal with two-nearby square pulses and a smaller additional one at the beginning.
+    n_baseline  = 64000
+    base_window = 1024
+    thr         = 5
+    wf          = np.zeros(n_baseline)
+    start1      = draw(integers(2*base_window   , n_baseline // 2))
+    length      = draw(integers(n_baseline // 50, n_baseline // 4))
+    stop1       = start1 + length
+    start2      = stop1 + 10
+    stop2       = start2 + length
+
+    heights = draw(lists(integers(thr+1, 50), min_size=3, max_size=3))
+    wf[start1:stop1] = np.full(length, heights[0])
+    wf[start2:stop2] = np.full(length, heights[1])
+    wf[1500  : 1520] = np.full(20    , heights[2]) # Extra pulse to check discard_width
+
+    return wf, thr, length, start1, start2
+
+
+@composite
+def trigger_time_and_validity(draw):
+    # Set of trigger times an validity flags, to test trigger coincidence
+    size  = draw (integers(min_value=1, max_value=100))
+    flag  = lists(integers(min_value=0, max_value=1  ), min_size=4, max_size=4)
+
+    i     = draw(lists(integers(min_value=0), min_size=size, max_size=size))
+    flags = draw(lists(flag                 , min_size=size, max_size=size))
+    return (i, flags)
+
+
+@given(double_square_wf=double_square_wf(), evt_id=integers(1, 10), pmt_id=integers(1, 10))
+def test_get_trigger_candidates_single(double_square_wf, evt_id, pmt_id):
+    # Given a square pulse check that is triggering where the pulse starts and values match the expected ones.
+    channel_config = channel_conf()
+    trigger_config = trigger_conf()
+
+    wf     = double_square_wf[0]
+    thr    = double_square_wf[1]
+    length = double_square_wf[2]
+    start1 = double_square_wf[3]
+    start2 = double_square_wf[4]
+
+    peak   = wf[start1:start2+length+channel_config['pulse_valid_ext']]
+    t_true = start1
+    q_true = peak[:-2].sum()
+    val_q  = in_range(q_true
+                     , channel_config['q_min'], channel_config['q_max']
+                     , left_closed=False, right_closed=False)
+    w_true = len(peak)
+    val_w  = in_range(w_true
+                     , channel_config['time_min'], channel_config['time_max']
+                     , left_closed=True , right_closed=False)
+    h_true = peak.max()
+    val_h  = h_true < channel_config['amp_max']
+    val_p  = True
+    val_t  = all([val_q, val_w, val_h, val_p])
+
+    triggers = get_trigger_candidates(wf, evt_id, pmt_id
+                                     ,channel_config, trigger_config
+                                     ,np.zeros(wf.shape), wf.max())
+
+    elements = [evt_id    , pmt_id
+               ,t_true    , q_true  , w_true, h_true
+               ,val_q     , val_w   , val_h , val_p , val_t
+               ,wf[start1-1], wf.max()
+               ,-1        , -1      , -1]
+
+    assert len(triggers) == 1
+    assert elements      == list(triggers[0])
+
+
+@given(double_square_wf=double_square_wf(), evt_id=integers(1, 10), pmt_id=integers(1, 10), pulse_valid=integers(0, 9))
+def test_get_trigger_candidates_double(double_square_wf, evt_id, pmt_id, pulse_valid):
+    # Check that the square pulses are identified as separated when pulse valid ext is lower than their separation.
+    channel_config = channel_conf(pulse_valid_ext=pulse_valid)
+    trigger_config = trigger_conf()
+
+    wf     = double_square_wf[0]
+    thr    = double_square_wf[1]
+    length = double_square_wf[2]
+    start1 = double_square_wf[3]
+    start2 = double_square_wf[4]
+
+    triggers = get_trigger_candidates(wf, evt_id, pmt_id
+                                     ,channel_config, trigger_config
+                                     ,np.zeros(wf.shape), wf.max())
+    assert len(triggers) == 2
+    for j, start in enumerate([start1, start2]):
+        peak   = wf[start:start+length+pulse_valid]
+        t_true = start
+        q_true = peak[:-2].sum()
+        val_q  = in_range(q_true
+                         , channel_config['q_min'], channel_config['q_max']
+                         , left_closed=False, right_closed=False)
+        w_true = len(peak)
+        val_w  = in_range(w_true
+                         , channel_config['time_min'], channel_config['time_max']
+                         , left_closed=True , right_closed=False)
+        h_true = peak.max()
+        val_h  = h_true < channel_config['amp_max']
+        val_p  = True
+        val_t  = all([val_q, val_w, val_h, val_p])
+
+        elements = [evt_id    , pmt_id
+                   ,t_true    , q_true  , w_true, h_true
+                   ,val_q     , val_w   , val_h , val_p , val_t
+                   ,wf[start1-1], wf.max()
+                   ,-1        , -1      , -1]
+
+        assert elements == list(triggers[j])
+
+
+@given(double_square_wf=double_square_wf(), evt_id=integers(1, 10), pmt_id=integers(1, 10), discard=integers(0, 30))
+def test_get_trigger_candidates_discard(double_square_wf, evt_id, pmt_id, discard):
+    # Check that the small pulse is kept when discard_width is low enough
+    channel_config = channel_conf()
+    trigger_config = trigger_conf(discard_width=discard)
+
+    wf     = double_square_wf[0]
+    thr    = double_square_wf[1]
+    length = double_square_wf[2]
+    start1 = double_square_wf[3]
+    start2 = double_square_wf[4]
+
+    triggers = get_trigger_candidates(wf, evt_id, pmt_id
+                                     ,channel_config, trigger_config
+                                     ,np.zeros(wf.shape), wf.max())
+    assert len(triggers) == 2
+    for j in range(2):
+        if j == 0:
+            peak   = wf[1500:1520+channel_config['pulse_valid_ext']]
+            t_true = 1500
+        else:
+            peak   = wf[start1:start2+length+channel_config['pulse_valid_ext']]
+            t_true = start1
+        q_true = peak[:-2].sum()
+        val_q  = in_range(q_true
+                         , channel_config['q_min'], channel_config['q_max']
+                         , left_closed=False, right_closed=False)
+        w_true = len(peak)
+        val_w  = in_range(w_true
+                         , channel_config['time_min'], channel_config['time_max']
+                         , left_closed=True , right_closed=False)
+        h_true = peak.max()
+        val_h  = h_true < channel_config['amp_max']
+        val_p  = True
+        val_t  = all([val_q, val_w, val_h, val_p])
+
+        elements = [evt_id    , pmt_id
+                   ,t_true    , q_true  , w_true, h_true
+                   ,val_q     , val_w   , val_h , val_p , val_t
+                   ,wf[start1-1], wf.max()
+                   ,-1        , -1      , -1]
+
+        assert elements == list(triggers[j])
+
+
+@given(double_square_wf=double_square_wf(), evt_id=integers(1, 10), pmt_id=integers(1, 10))
+def test_get_trigger_candidates_multipeak(double_square_wf, evt_id, pmt_id):
+    # Check that the multipeak protection discards the trigger if conditions are met.
+    channel_config = channel_conf(pulse_valid_ext=0)
+    trigger_config = trigger_conf()
+
+    wf     = double_square_wf[0]
+    thr    = double_square_wf[1]
+    length = double_square_wf[2]
+    start1 = double_square_wf[3]
+    start2 = double_square_wf[4]
+
+    triggers = get_trigger_candidates(wf, evt_id, pmt_id
+                                     ,channel_config, trigger_config
+                                     ,np.zeros(wf.shape), wf.max())
+    for t in triggers:
+        assert t.valid_all  == all(t[6:10])
+        assert t.valid_peak == True
+
+
+    ### Second peak characteristics
+    peak   = wf[start2:start2+length+channel_config['pulse_valid_ext']]
+    t_true = start2
+    q_true = peak[:-2].sum()
+    val_q  = in_range(q_true
+                     , channel_config['q_min'], channel_config['q_max']
+                     , left_closed=False, right_closed=False)
+    w_true = len(peak)
+    val_w  = in_range(w_true
+                     , channel_config['time_min'], channel_config['time_max']
+                     , left_closed=True , right_closed=False)
+    h_true = peak.max()
+    val_h  = h_true < channel_config['amp_max']
+    val_p  = True
+    val_t  = all([val_q, val_w, val_h, val_p])
+
+    # Multipeak condition metretrieve_trigger_information
+    trigger_config['multipeak'] = {'q_min':q_true//2, 'time_min':w_true//2, 'time_after':length+10}
+    triggers = get_trigger_candidates(wf, evt_id, pmt_id
+                                     ,channel_config, trigger_config
+                                     ,np.zeros(wf.shape), wf.max())
+
+    assert triggers[0].valid_peak == False
+    assert triggers[1].valid_peak == True
+    for t in triggers:
+        assert t.valid_all == all(t[6:10])
+
+    # Multipeak condition no longer met
+    trigger_config['multipeak'] = {'q_min':q_true*2, 'time_min':w_true*2, 'time_after':length+10}
+    triggers = get_trigger_candidates(wf, evt_id, pmt_id
+                                     ,channel_config, trigger_config
+                                     ,np.zeros(wf.shape), wf.max())
+
+    for t in triggers:
+        assert t.valid_all == all(t[6:10])
+        assert t.valid_peak == True
+
+
+@given(double_square_wf=double_square_wf(), evt_id=integers(0, 3))
+def test_retrieve_trigger_information(double_square_wf, evt_id):
+    channel_config = channel_conf()
+    trigger_config = trigger_conf()
+
+    retriever = retrieve_trigger_information({0:channel_config, 1:channel_config}, trigger_config)
+
+    wf     = double_square_wf[0]
+    wfs    = np.array([i*wf for i in range(1, 4)])
+    thr    = double_square_wf[1]
+    length = double_square_wf[2]
+    start1 = double_square_wf[3]
+    start2 = double_square_wf[4]
+
+    triggers  = retriever(-wfs, wfs[:2], np.zeros(wfs.shape), evt_id)
+
+    assert len(triggers) == 2
+    for j in range(2):
+        peak   = wfs[j][start1:start2+length+channel_config['pulse_valid_ext']]
+        t_true = start1
+        q_true = peak[:-2].sum()
+        val_q  = in_range(q_true
+                         , channel_config['q_min'], channel_config['q_max']
+                         , left_closed=False, right_closed=False)
+        w_true = len(peak)
+        val_w  = in_range(w_true
+                         , channel_config['time_min'], channel_config['time_max']
+                         , left_closed=True , right_closed=False)
+        h_true = peak.max()
+        val_h  = h_true < channel_config['amp_max']
+        val_p  = True
+        val_t  = all([val_q, val_w, val_h, val_p])
+
+        elements = [evt_id    , j
+                   ,t_true    , q_true  , w_true, h_true
+                   ,val_q     , val_w   , val_h , val_p , val_t
+                   ,wf[start1-1], -wfs[j].max()
+                   ,-1        , -1      , -1]
+
+        assert np.allclose(elements, list(triggers[j]))
+
+
+@given(trigger_time_and_validity(), integers())
+def test_check_trigger_coincidence(trigger_info, coincidence_window):
+    trigger_time = trigger_info[0]
+    validity     = trigger_info[1]
+    trigg        = np.zeros((len(validity), 16), dtype='object')
+
+    # Initiate sim. trigger values
+    trigg[:, 1   ] = list(range(len(validity)))
+    trigg[:, 2   ] = trigger_time
+    trigg[:, 6:11] = True#validity
+    trigg[:, -3: ] = np.full((len(validity), 3), -1)
+
+    order_idx  = np.lexsort((trigg[:,1], trigger_time))
+    old_order  = np.argsort(order_idx)
+
+    ordered_times = trigg[:, 2   ][order_idx]
+    ordered_val   = trigg[:, 6:11][order_idx]
+    ordered_ids   = trigg[:, 1   ][order_idx]
+    trigger_back  = trigg[order_idx].copy()
+
+    for i, time1 in enumerate(ordered_times):
+        if not np.all(ordered_val[i]): continue
+        trigger_back[i, -3] = 1
+        nearest = False
+        for val2, id2, time2 in zip(ordered_val[i+1:], ordered_ids[i+1:], ordered_times[i+1:]):
+            if not np.all(val2): continue
+            if not nearest:
+                trigger_back[i, -2] = time2-time1
+                trigger_back[i, -1] = id2
+                nearest = True
+            if (time2-time1)<coincidence_window:
+                trigger_back[i, -3] += 1
+            else: break
+
+    checks_coinc = check_trigger_coincidence(coincidence_window)
+    new_t        = checks_coinc([TriggerInfo(*t) for t in trigg])
+
+    for t_old, t_new in zip(trigger_back[old_order], new_t):
+        assert list(t_old) == list(t_new)
diff --git a/invisible_cities/sierpe/blr.pxd b/invisible_cities/sierpe/blr.pxd
index ad6fbd15c5..76710c3546 100644
--- a/invisible_cities/sierpe/blr.pxd
+++ b/invisible_cities/sierpe/blr.pxd
@@ -32,4 +32,10 @@ cpdef deconvolve_signal(double [:] signal_daq,
                         double     coeff_clean     = *,
                         double     coeff_blr       = *,
                         double     thr_trigger     = *,
-                        int accum_discharge_length = *)
\ No newline at end of file
+                        int accum_discharge_length = *)
+
+cpdef deconvolve_signal_fpga(short [:] signal_daq
+                            ,double    coeff_clean = *
+                            ,double    coeff_blr   = *
+                            ,double    thr_trigger = *
+                            ,size_t    base_window = *)
diff --git a/invisible_cities/sierpe/blr.pyx b/invisible_cities/sierpe/blr.pyx
index 3db7863127..279da44f94 100644
--- a/invisible_cities/sierpe/blr.pyx
+++ b/invisible_cities/sierpe/blr.pyx
@@ -1,3 +1,4 @@
+#cython: language_level=3
 import  numpy as np
 cimport numpy as np
 from scipy import signal as SGN
@@ -59,7 +60,6 @@ cpdef deconvolve_signal(double [:] signal_daq,
 
         if (signal_daq[k] < trigger_line) and (acum[k-1] < thr_acum):
             # discharge accumulator
-
             if acum[k-1] > 1:
                 acum[k] = acum[k-1] * (1 - coef)
                 if j < accum_discharge_length - 1:
@@ -70,4 +70,138 @@ cpdef deconvolve_signal(double [:] signal_daq,
                 acum[k] = 0
                 j = 0
     # return recovered signal
-    return np.asarray(signal_r)
\ No newline at end of file
+    return np.asarray(signal_r)
+
+
+cpdef deconvolve_signal_fpga( short [:] signal_daq
+                            , double    coeff_clean = 2.905447E-06
+                            , double    coeff_blr   = 1.632411E-03
+                            , double    thr_trigger = 5
+                            , size_t    base_window = 1024        ):
+
+    """
+    Simulate the deconvolution process as in the daq, differences compared to
+    usual offline deconvolution:
+      - Baseline is calculated as a moving average of 1024 counts (FPGA).
+      - Flips result between raw and deconvolved signal when outside of both
+        signal and discharge regions.
+      - Discharge made at fixed value (0.995)
+      - Result is truncated to integers.
+      - ADC threshold acting as absolute threshold.
+
+    Parameters
+    ----------
+    signal_daq  : short array
+         PMT raw waveform
+    coeff_clean : double
+         Characteristic parameter of the high pass filter
+    coeff_blr   : double
+         Characteristic parameter of BLR
+    thr_trigger : double
+         Threshold in ADCs to activate BLR
+    base_window : size
+         Moving average window for baseline calculation
+
+    Returns
+    ----------
+    Tuple of arrays:
+    - Deconvolved waveform (int16  [:])
+    - Baseline             (int16  [:])
+    """
+    cdef double coef = coeff_blr
+    cdef double thr_acum = thr_trigger / coef
+    cdef int len_signal_daq = len(signal_daq)
+
+    cdef double [:] signal_r = np.zeros(len_signal_daq, dtype=np.double)
+    cdef double [:] signal_f = np.zeros(len_signal_daq, dtype=np.double)
+    cdef double [:] acum     = np.zeros(len_signal_daq, dtype=np.double)
+
+    cdef int j
+
+    # compute noise
+    cdef double noise =  0
+    cdef int nn = 400 # fixed at 10 mus
+
+    for j in range(nn):
+        noise += signal_daq[j] * signal_daq[j]
+    noise /= nn
+    cdef double noise_rms = np.sqrt(noise)
+
+    # trigger line
+    cdef double trigger_line = thr_trigger #* noise_rms
+    # cleaning signal
+    cdef double [:]  b_cf
+    cdef double [:]  a_cf
+
+    ### Baseline related variables
+    cdef double [:]  top   = np.zeros(len_signal_daq, dtype=np.double)
+    cdef short  [:]  aux   = np.zeros(len_signal_daq, dtype=np.int16)
+    cdef long ped          = np.sum  (signal_daq[0:base_window], dtype=np.int32)
+    cdef size_t delay      = 2 # Delay in the FPGA in the bin used for baseline substraction
+    cdef unsigned int iaux = base_window
+
+    top[0:base_window] = ped/base_window
+    aux[0:base_window] = signal_daq[0:base_window]
+
+    b_cf, a_cf = SGN.butter(1, coeff_clean, 'high', analog=False);
+    g, a1 = b_cf[0], -a_cf[-1]
+
+    ### Initiate filt
+    filt_hr   = 0
+    #1st
+    filt_x    = g  * -(signal_daq[0] - top[0])
+    filt_a1hr = a1 * filt_hr
+
+    filt_h    = filt_x + filt_a1hr
+    signal_f[0] = filt_h - filt_hr
+    filt_hr     = filt_h
+
+    #2nd
+    filt_x      = g  * -(signal_daq[1] - top[0])
+    filt_a1hr   = a1 * filt_hr
+
+    filt_h      = filt_x + filt_a1hr
+    signal_f[1] = filt_h - filt_hr
+    filt_hr     = filt_h
+
+
+    cdef int k
+
+    #Initiate BLR
+    for k in range(0, delay):
+        signal_r[k] = signal_daq[k]
+
+    for k in range(2, len_signal_daq):
+        # always update signal and accumulator
+        current     = signal_daq[k]
+        baseline    = top[k-delay]
+
+        ### High-pass filter
+        filt_x      = g  * -(current - baseline)
+        filt_a1hr   = a1 * filt_hr
+
+        filt_h      = filt_x + filt_a1hr
+        signal_f[k] = filt_h - filt_hr
+
+        ### BLR restoration
+        blr_val     = (signal_f[k] + signal_f[k]*(coef / 2) +
+                       coef * acum[k-1]) + baseline
+
+        acum[k] = acum[k-1] + signal_f[k]
+
+        if (signal_f[k] < trigger_line) and (acum[k-1] < thr_acum):
+            # discharge accumulator
+            if acum[k-1] > 1:
+                acum[k] = acum[k-1] * .995 #Fixed discharge in FPGA code
+            else:
+                acum [k] = 0
+                blr_val = current # When outside of BLR/discharge, flip to raw signal
+                if k>=base_window:
+                    ped       = ped + current - aux[iaux-base_window]
+                    aux[iaux] = current
+                    iaux     += 1
+
+        signal_r[k] = blr_val
+        top[k]      = ped/base_window
+
+    return np.asarray(signal_r, dtype='int16'), np.asarray(top, dtype='int16')
diff --git a/invisible_cities/sierpe/blr_test.py b/invisible_cities/sierpe/blr_test.py
index 7e5b01a055..5bc6c9351e 100644
--- a/invisible_cities/sierpe/blr_test.py
+++ b/invisible_cities/sierpe/blr_test.py
@@ -2,19 +2,25 @@
 
 import numpy  as np
 import tables as tb
+import pandas as pd
 
-from pytest import fixture
-from pytest import mark
-from flaky  import flaky
+from pytest                    import fixture
+from pytest                    import mark
+from flaky                     import flaky
+from hypothesis                import given
+from hypothesis                import settings
+from hypothesis.strategies     import integers
 
 from .. reco import calib_sensors_functions as csf
 from .       import blr
 
 
-deconv_params = namedtuple("deconv_params",
-                           "coeff_clean coeff_blr "
-                           "thr_trigger accum_discharge_length")
-
+deconv_params      = namedtuple("deconv_params",
+                                "coeff_clean coeff_blr "
+                                "thr_trigger accum_discharge_length")
+deconv_fpga_params = namedtuple("deconv_params",
+                                "coeff_clean coeff_blr "
+                                "thr_trigger base_window")
 
 @fixture(scope="session")
 def sin_wf_params():
@@ -54,6 +60,20 @@ def ad_hoc_blr_signals(example_blr_wfs_filename):
         return rwf, blrwf, attrs.n_baseline, params
 
 
+@fixture(scope="session")
+def ad_hoc_blr_fpga_signals(example_blr_fpga_wfs_filename):
+    with tb.open_file(example_blr_fpga_wfs_filename) as file:
+        rwf    = file.root.RD.pmtrwf[:]
+        blrwf  = file.root.BLR[:]
+        base   = file.root.Baseline[:]
+        attrs  = file.root.BLR.attrs
+        params = deconv_fpga_params(attrs.coeff_c    ,
+                                    attrs.coeff_blr  ,
+                                    attrs.thr_trigger,
+                                    attrs.base_window)
+        return rwf, blrwf, base, params
+
+
 def test_deconvolve_signal_positive_integral(sin_wf, sin_wf_params):
     # The RWF should have null integral because contains roughly
     # the same number of positive and negative samples.
@@ -90,7 +110,7 @@ def test_deconvolve_signal_ad_hoc_signals(ad_hoc_blr_signals):
 
     rwf         = all_rwfs        [evt_no, pmt_no]
     true_blr_wf = all_true_blr_wfs[evt_no, pmt_no]
-    
+
     cwf         = np.mean(rwf[:n_baseline]) - rwf
     blr_wf = blr.deconvolve_signal(cwf,
                                    coeff_clean            = params.coeff_clean[pmt_no],
@@ -144,3 +164,62 @@ def test_deconv_pmt_ad_hoc_signals_dead_sensors(ad_hoc_blr_signals):
                                   rep_thr              , rep_acc             )))
 
     np.allclose(blr_wfs, evt_true_blr_wfs[pmt_active])
+
+
+@settings(max_examples=1, deadline=1000)
+@given   (evt_no=integers(0, 9))
+@mark.slow
+def test_deconv_pmt_fpga_ad_hoc_signals_all(evt_no, ad_hoc_blr_fpga_signals):
+    all_rwfs, all_true_blr_wfs, all_true_baseline, params = ad_hoc_blr_fpga_signals
+
+    evt_rwfs          = all_rwfs         [evt_no]
+    evt_true_blr_wfs  = all_true_blr_wfs [evt_no]
+    evt_true_baseline = all_true_baseline[evt_no]
+
+    rep_thr  = np.repeat(params.thr_trigger, evt_rwfs.shape[0])
+    rep_base = np.repeat(params.base_window, evt_rwfs.shape[0])
+    blr_wfs, baseline = zip(*map(blr.deconvolve_signal_fpga, evt_rwfs        ,
+                                 params.coeff_clean        , params.coeff_blr,
+                                 rep_thr                   , rep_base        ))
+
+    assert np.allclose(blr_wfs , evt_true_blr_wfs )
+    assert np.allclose(baseline, evt_true_baseline)
+
+
+@settings(max_examples=1, deadline=1000)
+@given   (evt_no=integers(0, 9))
+@mark.slow
+def test_deconv_pmt_fpga_ad_hoc_signals_dead_sensors(evt_no, ad_hoc_blr_fpga_signals):
+    all_rwfs, all_true_blr_wfs, all_true_baseline, params = ad_hoc_blr_fpga_signals
+
+    n_evts, n_pmts, _ = all_rwfs.shape
+    pmt_active        = np.arange(n_pmts)
+    n_alive           = np.random.randint(1, n_pmts - 1)
+    pmt_active        = np.random.choice(pmt_active, size=n_alive, replace=False)
+
+    evt_rwfs          = all_rwfs         [evt_no]
+    evt_true_blr_wfs  = all_true_blr_wfs [evt_no]
+    evt_true_baseline = all_true_baseline[evt_no]
+
+    rep_thr  = np.repeat(params.thr_trigger, len(pmt_active))
+    rep_base = np.repeat(params.base_window, len(pmt_active))
+    blr_wfs, baseline = zip(*map(blr.deconvolve_signal_fpga    , evt_rwfs        [pmt_active],
+                                 params.coeff_clean[pmt_active], params.coeff_blr[pmt_active],
+                                 rep_thr                       , rep_base                    ))
+
+    assert np.allclose(blr_wfs , evt_true_blr_wfs [pmt_active])
+    assert np.allclose(baseline, evt_true_baseline[pmt_active])
+
+
+@given(thr=integers(2, 100), window=integers(10, 5000))
+def test_moving_average_baseline(thr, window):
+    samples = 10000
+    np.random.seed(1234)
+    wvf     = np.random.randint(0 , thr-1, samples, dtype=np.int16)
+
+    _, baseline = blr.deconvolve_signal_fpga(wvf, 1e-6, 1e-3, thr, window)
+
+    rolling = pd.Series(wvf).rolling(window=window).mean().values.astype(np.int16)
+
+    assert np.allclose(rolling[window:], baseline[window:])
+    assert np.allclose(np.repeat(wvf[:window].mean().astype(np.int16), window), baseline[:window])