diff --git a/libraries/base/common/hdllib.cfg b/libraries/base/common/hdllib.cfg
index 75b70901a1163728a9d1d7f6d3f551db58b6384d..e95d804fc699c49196b5e95764f559a01854ca91 100644
--- a/libraries/base/common/hdllib.cfg
+++ b/libraries/base/common/hdllib.cfg
@@ -8,9 +8,9 @@ synth_files =
     $RADIOHDL_WORK/libraries/base/common/src/vhdl/common_pkg.vhd
     src/vhdl/common_str_pkg.vhd
     src/vhdl/common_mem_pkg.vhd
-    src/vhdl/common_math_pkg.vhd
     src/vhdl/common_field_pkg.vhd
     src/vhdl/common_lfsr_sequences_pkg.vhd
+    src/vhdl/common_math_pkg.vhd
     src/vhdl/common_interface_layers_pkg.vhd
     src/vhdl/common_network_layers_pkg.vhd
     src/vhdl/common_network_total_header_pkg.vhd
diff --git a/libraries/base/common/src/vhdl/common_math_pkg.vhd b/libraries/base/common/src/vhdl/common_math_pkg.vhd
index 0cbb134183563d9fc256e750614fca2a017625a0..1fd31bbd370afc179e24e0247d8a8f091e12df39 100644
--- a/libraries/base/common/src/vhdl/common_math_pkg.vhd
+++ b/libraries/base/common/src/vhdl/common_math_pkg.vhd
@@ -36,7 +36,7 @@ LIBRARY IEEE;
 USE IEEE.STD_LOGIC_1164.ALL;
 USE IEEE.MATH_REAL.ALL;
 USE work.common_pkg.ALL;
-
+USE work.common_lfsr_sequences_pkg.ALL;
 
 PACKAGE common_math_pkg IS
 
@@ -79,7 +79,10 @@ PACKAGE common_math_pkg IS
   -- To create an FFT input phasor with frequency in the middle of a channel use FREQ = ch.
   FUNCTION common_math_create_look_up_table_phasor(N, W : POSITIVE; AMPL, FREQ, PHI : REAL) RETURN t_nat_integer_arr;
   FUNCTION common_math_create_look_up_table_phasor(N, W : POSITIVE; AMPL, FREQ, PHI : REAL) RETURN t_slv_32_arr;  -- range 0 TO N-1
-  
+ 
+  FUNCTION common_math_create_random_arr(N, W : POSITIVE; seed : NATURAL) RETURN t_integer_arr;
+
+ 
 END common_math_pkg;
 
 
@@ -178,6 +181,18 @@ PACKAGE BODY common_math_pkg IS
     END LOOP;
     RETURN v_exp_arr;
   END;
-  
+   
+  FUNCTION common_math_create_random_arr(N, W : POSITIVE; seed : NATURAL) RETURN t_integer_arr IS
+    VARIABLE v_rand_arr : t_integer_arr(0 TO N-1);
+    VARIABLE v_random : STD_LOGIC_VECTOR(W-1 DOWNTO 0) := TO_UVEC(seed, W);
+  BEGIN
+    FOR I IN 0 TO N-1 LOOP
+      v_random := func_common_random(v_random);
+      v_rand_arr(I) := TO_SINT(v_random);
+    END LOOP;
+    RETURN v_rand_arr;
+  END;
+
+
 END common_math_pkg;
 
diff --git a/libraries/dsp/st/hdllib.cfg b/libraries/dsp/st/hdllib.cfg
index 3617a5976a7da4b76d245c3cbcd1c1cfdcfb6d71..faba09d3db7d4382948bb971bffd4a62e111bcd4 100644
--- a/libraries/dsp/st/hdllib.cfg
+++ b/libraries/dsp/st/hdllib.cfg
@@ -9,22 +9,29 @@ synth_files =
     src/vhdl/st_ctrl.vhd 
     src/vhdl/st_calc.vhd 
     src/vhdl/st_sst.vhd 
+    src/vhdl/st_xsq.vhd 
+    src/vhdl/st_xsq_arr.vhd 
 #    src/vhdl/st_top.vhd 
     src/vhdl/st_histogram.vhd
     src/vhdl/st_histogram_reg.vhd
     src/vhdl/mms_st_histogram.vhd
     src/vhdl/st_histogram_8_april.vhd
+
+    tb/vhdl/tb_st_pkg.vhd 
  
 test_bench_files = 
     tb/vhdl/tb_st_acc.vhd 
     tb/vhdl/tb_st_calc.vhd 
     tb/vhdl/tb_mmf_st_sst.vhd
+    tb/vhdl/tb_st_xsq.vhd
+    tb/vhdl/tb_tb_st_xsq.vhd
     tb/vhdl/tb_st_histogram.vhd
     tb/vhdl/tb_mms_st_histogram.vhd
     tb/vhdl/tb_tb_st_histogram.vhd
 
 regression_test_vhdl = 
     tb/vhdl/tb_st_acc.vhd 
+    tb/vhdl/tb_tb_st_xsq.vhd
     #tb/vhdl/tb_st_calc.vhd   -- tb is not self checking yet
 
 
diff --git a/libraries/dsp/st/src/vhdl/st_xsq.vhd b/libraries/dsp/st/src/vhdl/st_xsq.vhd
new file mode 100644
index 0000000000000000000000000000000000000000..997e13c1b253434c103dca1bc505447b0bfa37ab
--- /dev/null
+++ b/libraries/dsp/st/src/vhdl/st_xsq.vhd
@@ -0,0 +1,255 @@
+-------------------------------------------------------------------------------
+--
+-- Copyright 2021
+-- ASTRON (Netherlands Institute for Radio Astronomy) <http://www.astron.nl/>
+-- P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
+--
+-- Licensed under the Apache License, Version 2.0 (the "License");
+-- you may not use this file except in compliance with the License.
+-- You may obtain a copy of the License at
+--
+--     http://www.apache.org/licenses/LICENSE-2.0
+--
+-- Unless required by applicable law or agreed to in writing, software
+-- distributed under the License is distributed on an "AS IS" BASIS,
+-- WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+-- See the License for the specific language governing permissions and
+-- limitations under the License.
+--
+-------------------------------------------------------------------------------
+-- Author : R. vd Walle
+-- Purpose:
+--   Store the statistics of the complex input streams in_a and in_b with
+--   blocks of g_nof_crosslets * g_nof_signal_inputs**2 multiplexed subbands into a MM register.
+-- Description:                                                          
+--
+--   After each sync the MM register gets updated with the statistics
+--   of the previous sync interval. The length of the sync interval determines
+--   the nof accumlations per statistic, hence the integration time. See st_calc
+--   for more details.
+-- Remarks:
+-- . The in_sync is assumed to be a pulse an interpreted directly.
+-- . The MM register is single page RAM to save memory resources. Therefore
+--   just after the sync its contents is undefined when it gets written, but
+--   after that its contents remains stable for the rest of the sync interval.
+--   Therefore it is not necessary to use a dual page register that swaps at
+--   the sync. 
+-- . The minimum c_nof_statistics = 8. Lower values lead to simulation errors. This is
+--   due to the read latency of 2 of the accumulation memory in the st_calc entity.
+-- . More detail can be found in:
+--   https://support.astron.nl/confluence/display/L2M/L5+SDPFW+Design+Document%3A+Subband+Correlator
+-------------------------------------------------------------------------------
+
+LIBRARY IEEE, common_lib, mm_lib, technology_lib, dp_lib;
+USE IEEE.std_logic_1164.ALL;
+USE common_lib.common_pkg.ALL;
+USE common_lib.common_mem_pkg.ALL;
+USE common_lib.common_field_pkg.ALL;
+USE dp_lib.dp_stream_pkg.ALL;
+USE technology_lib.technology_select_pkg.ALL;
+ENTITY st_xsq IS
+  GENERIC (
+    g_nof_signal_inputs : NATURAL := 2;
+    g_nof_crosslets     : NATURAL := 1;
+    g_in_data_w         : NATURAL := 18;    -- width of the data to be accumulated
+    g_stat_data_w       : NATURAL := 54;    -- statistics accumulator width
+    g_stat_data_sz      : NATURAL := 2      -- statistics word width >= statistics accumulator width and fit in a power of 2 multiple 32b MM words
+  );                
+  PORT (            
+    mm_rst : IN  STD_LOGIC;
+    mm_clk : IN  STD_LOGIC;
+    dp_rst : IN  STD_LOGIC;
+    dp_clk : IN  STD_LOGIC;
+                    
+    -- Streaming    
+    in_a   : IN  t_dp_sosi;   -- Complex input data
+    in_b   : IN  t_dp_sosi;   -- Complex input data
+    
+    -- Memory Mapped
+    ram_st_xsq_mosi : IN  t_mem_mosi := c_mem_mosi_rst;  
+    ram_st_xsq_miso : OUT t_mem_miso := c_mem_miso_rst
+  );
+END st_xsq;
+
+ARCHITECTURE str OF st_xsq IS
+
+  CONSTANT c_xsq               : NATURAL := g_nof_signal_inputs * g_nof_signal_inputs;
+  CONSTANT c_nof_statistics    : NATURAL := g_nof_crosslets * c_xsq;
+  CONSTANT c_nof_stat_w        : NATURAL := ceil_log2(c_nof_statistics);
+  CONSTANT c_nof_word          : NATURAL := g_stat_data_sz*c_nof_statistics;
+  CONSTANT c_nof_word_w        : NATURAL := ceil_log2(c_nof_word);
+  CONSTANT c_stat_word_w       : NATURAL := g_stat_data_sz*c_word_w;
+  CONSTANT c_total_ram_addr_w  : NATURAL := ceil_log2(c_nof_complex) + c_nof_word_w;
+  
+  -- Statistics register
+  CONSTANT c_mm_ram   : t_c_mem := (latency  => 1,
+                                    adr_w    => c_nof_word_w,
+                                    dat_w    => c_word_w,
+                                    nof_dat  => c_nof_word,
+                                    init_sl  => '0');           -- MM side : sla_in, sla_out
+  CONSTANT c_stat_ram : t_c_mem := (latency  => 1,
+                                    adr_w    => c_nof_stat_w,
+                                    dat_w    => c_stat_word_w,
+                                    nof_dat  => c_nof_statistics,
+                                    init_sl  => '0');           -- ST side : stat_mosi
+  SIGNAL pipe_in_a    : t_dp_sosi;
+  SIGNAL pipe_in_b    : t_dp_sosi;  
+  
+  SIGNAL stat_data_re : STD_LOGIC_VECTOR(g_stat_data_w-1 DOWNTO 0);  
+  SIGNAL stat_data_im : STD_LOGIC_VECTOR(g_stat_data_w-1 DOWNTO 0);  
+  
+  SIGNAL wrdata_re    : STD_LOGIC_VECTOR(c_mem_data_w-1 DOWNTO 0);
+  SIGNAL wrdata_im    : STD_LOGIC_VECTOR(c_mem_data_w-1 DOWNTO 0);
+  
+  SIGNAL stat_mosi    : t_mem_mosi := c_mem_mosi_rst;
+
+  SIGNAL ram_st_xsq_mosi_arr      : t_mem_mosi_arr(c_nof_complex-1 DOWNTO 0) := (OTHERS => c_mem_mosi_rst);
+  SIGNAL ram_st_xsq_miso_arr      : t_mem_miso_arr(c_nof_complex-1 DOWNTO 0) := (OTHERS => c_mem_miso_rst);
+  SIGNAL remapped_ram_st_xsq_mosi : t_mem_mosi;  
+    
+BEGIN
+  ---------------------------------------------------------------
+  -- pipeline inputs to increase latency with 1 in comparison to sync for st_calc 
+  ---------------------------------------------------------------
+  u_dp_pipeline_a : ENTITY dp_lib.dp_pipeline
+  GENERIC MAP (
+    g_pipeline => 1
+  )
+  PORT MAP (
+    rst          => dp_rst,
+    clk          => dp_clk,
+    -- ST sink
+    snk_in       => in_a,
+    -- ST source
+    src_out      => pipe_in_a
+  );
+
+  u_dp_pipeline_b : ENTITY dp_lib.dp_pipeline
+  GENERIC MAP (
+    g_pipeline => 1
+  )
+  PORT MAP (
+    rst          => dp_rst,
+    clk          => dp_clk,
+    -- ST sink
+    snk_in       => in_b,
+    -- ST source
+    src_out      => pipe_in_b
+  );
+
+  ---------------------------------------------------------------
+  -- st_calc 
+  ---------------------------------------------------------------
+  st_calc : ENTITY work.st_calc 
+  GENERIC MAP (
+    g_nof_mux      => 1,
+    g_nof_stat     => c_nof_statistics,
+    g_in_dat_w     => g_in_data_w,
+    g_out_dat_w    => g_stat_data_w,
+    g_out_adr_w    => c_nof_stat_w,
+    g_complex      => TRUE
+  )
+  PORT MAP (
+    rst            => dp_rst,
+    clk            => dp_clk,
+    in_ar          => pipe_in_a.re(g_in_data_w-1 DOWNTO 0),
+    in_ai          => pipe_in_a.im(g_in_data_w-1 DOWNTO 0),
+    in_br          => pipe_in_b.re(g_in_data_w-1 DOWNTO 0),
+    in_bi          => pipe_in_b.im(g_in_data_w-1 DOWNTO 0),
+    in_val         => pipe_in_a.valid,
+    in_sync        => in_a.sync,
+    out_adr        => stat_mosi.address(c_stat_ram.adr_w-1 DOWNTO 0),
+    out_re         => stat_data_re,
+    out_im         => stat_data_im,
+    out_val        => stat_mosi.wr,
+    out_val_m      => OPEN
+  );
+  
+  wrdata_re <= RESIZE_MEM_UDATA(stat_data_re);
+  wrdata_im <= RESIZE_MEM_UDATA(stat_data_im);
+  
+  ---------------------------------------------------------------
+  -- COMBINE MEMORY MAPPED INTERFACES
+  ---------------------------------------------------------------
+  -- Translate incoming MM interface [crosslets][in A][in B][complex][word] to
+  -- [complex][crosslets][in A][in B][word]
+  p_remap : PROCESS(ram_st_xsq_mosi)
+  BEGIN
+    remapped_ram_st_xsq_mosi <= ram_st_xsq_mosi;
+    remapped_ram_st_xsq_mosi.address(c_total_ram_addr_w -1 DOWNTO c_nof_word_w) <= ram_st_xsq_mosi.address(ceil_log2(c_nof_complex)+ceil_log2(g_stat_data_sz)-1 DOWNTO ceil_log2(g_stat_data_sz));
+    remapped_ram_st_xsq_mosi.address(c_nof_word_w -1 DOWNTO 0) <= ram_st_xsq_mosi.address(c_total_ram_addr_w -1 DOWNTO ceil_log2(c_nof_complex)+ceil_log2(g_stat_data_sz)) & ram_st_xsq_mosi.address(ceil_log2(g_stat_data_sz)-1 DOWNTO 0);
+  END PROCESS;
+
+  -- Combine the internal array of mm interfaces for both real
+  -- and imaginary part. 
+  u_mem_mux_select : entity common_lib.common_mem_mux
+  generic map (    
+    g_nof_mosi    => c_nof_complex,
+    g_mult_addr_w => c_nof_word_w
+  )
+  port map (
+    mosi     => remapped_ram_st_xsq_mosi,
+    miso     => ram_st_xsq_miso,
+    mosi_arr => ram_st_xsq_mosi_arr,
+    miso_arr => ram_st_xsq_miso_arr
+  );
+ 
+  -- ram for real values 
+  stat_reg_re : ENTITY common_lib.common_ram_crw_crw_ratio
+  GENERIC MAP (
+    g_ram_a      => c_mm_ram,
+    g_ram_b      => c_stat_ram,
+    g_init_file  => "UNUSED"
+  )
+  PORT MAP (
+    rst_a     => mm_rst,
+    clk_a     => mm_clk,
+    
+    rst_b     => dp_rst,
+    clk_b     => dp_clk,
+    
+    wr_en_a   => ram_st_xsq_mosi_arr(0).wr,  -- only for diagnostic purposes, typically statistics are read only
+    wr_dat_a  => ram_st_xsq_mosi_arr(0).wrdata(c_mm_ram.dat_w-1 DOWNTO 0),
+    adr_a     => ram_st_xsq_mosi_arr(0).address(c_mm_ram.adr_w-1 DOWNTO 0),
+    rd_en_a   => ram_st_xsq_mosi_arr(0).rd,
+    rd_dat_a  => ram_st_xsq_miso_arr(0).rddata(c_mm_ram.dat_w-1 DOWNTO 0),
+    rd_val_a  => ram_st_xsq_miso_arr(0).rdval,
+    
+    wr_en_b   => stat_mosi.wr,
+    wr_dat_b  => wrdata_re(c_stat_ram.dat_w-1 DOWNTO 0),
+    adr_b     => stat_mosi.address(c_stat_ram.adr_w-1 DOWNTO 0),
+    rd_en_b   => '0',
+    rd_dat_b  => OPEN,
+    rd_val_b  => OPEN
+  );                                   
+ 
+  -- ram for imaginary values
+  stat_reg_im : ENTITY common_lib.common_ram_crw_crw_ratio
+  GENERIC MAP (
+    g_ram_a      => c_mm_ram,
+    g_ram_b      => c_stat_ram,
+    g_init_file  => "UNUSED"
+  )
+  PORT MAP (
+    rst_a     => mm_rst,
+    clk_a     => mm_clk,
+    
+    rst_b     => dp_rst,
+    clk_b     => dp_clk,
+    
+    wr_en_a   => ram_st_xsq_mosi_arr(1).wr,  -- only for diagnostic purposes, typically statistics are read only
+    wr_dat_a  => ram_st_xsq_mosi_arr(1).wrdata(c_mm_ram.dat_w-1 DOWNTO 0),
+    adr_a     => ram_st_xsq_mosi_arr(1).address(c_mm_ram.adr_w-1 DOWNTO 0),
+    rd_en_a   => ram_st_xsq_mosi_arr(1).rd,
+    rd_dat_a  => ram_st_xsq_miso_arr(1).rddata(c_mm_ram.dat_w-1 DOWNTO 0),
+    rd_val_a  => ram_st_xsq_miso_arr(1).rdval,
+    
+    wr_en_b   => stat_mosi.wr,
+    wr_dat_b  => wrdata_im(c_stat_ram.dat_w-1 DOWNTO 0),
+    adr_b     => stat_mosi.address(c_stat_ram.adr_w-1 DOWNTO 0),
+    rd_en_b   => '0',
+    rd_dat_b  => OPEN,
+    rd_val_b  => OPEN
+  );     
+  
+END str;
diff --git a/libraries/dsp/st/src/vhdl/st_xsq_arr.vhd b/libraries/dsp/st/src/vhdl/st_xsq_arr.vhd
new file mode 100644
index 0000000000000000000000000000000000000000..b4445b2022a4ef884e3c5eeb16bf8d465ff683ac
--- /dev/null
+++ b/libraries/dsp/st/src/vhdl/st_xsq_arr.vhd
@@ -0,0 +1,114 @@
+-------------------------------------------------------------------------------
+--
+-- Copyright 2021
+-- ASTRON (Netherlands Institute for Radio Astronomy) <http://www.astron.nl/>
+-- P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
+--
+-- Licensed under the Apache License, Version 2.0 (the "License");
+-- you may not use this file except in compliance with the License.
+-- You may obtain a copy of the License at
+--
+--     http://www.apache.org/licenses/LICENSE-2.0
+--
+-- Unless required by applicable law or agreed to in writing, software
+-- distributed under the License is distributed on an "AS IS" BASIS,
+-- WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+-- See the License for the specific language governing permissions and
+-- limitations under the License.
+--
+-------------------------------------------------------------------------------
+-- Author : R. vd Walle
+-- Purpose:
+-- Multiple instances of st_xsq.vhd
+-- Description:                                                          
+--   . See st_xsq.vhd
+-- Remarks:
+-------------------------------------------------------------------------------
+
+LIBRARY IEEE, common_lib, mm_lib, technology_lib, dp_lib;
+USE IEEE.std_logic_1164.ALL;
+USE common_lib.common_pkg.ALL;
+USE common_lib.common_mem_pkg.ALL;
+USE common_lib.common_field_pkg.ALL;
+USE dp_lib.dp_stream_pkg.ALL;
+ENTITY st_xsq_arr IS
+  GENERIC (
+    g_nof_streams       : NATURAL := 1;
+    g_nof_crosslets     : NATURAL := 1;
+    g_nof_signal_inputs : NATURAL := 2;
+    g_in_data_w         : NATURAL := 18;    -- width of the data to be accumulated
+    g_stat_data_w       : NATURAL := 54;    -- statistics accumulator width
+    g_stat_data_sz      : NATURAL := 2      -- statistics word width >= statistics accumulator width and fit in a power of 2 multiple 32b MM words
+  );                
+  PORT (            
+    mm_rst : IN  STD_LOGIC;
+    mm_clk : IN  STD_LOGIC;
+    dp_rst : IN  STD_LOGIC;
+    dp_clk : IN  STD_LOGIC;
+                    
+    -- Streaming    
+    in_a_arr : IN  t_dp_sosi_arr(g_nof_streams-1 DOWNTO 0);   -- Complex input data
+    in_b_arr : IN  t_dp_sosi_arr(g_nof_streams-1 DOWNTO 0);   -- Complex input data
+    
+    -- Memory Mapped
+    ram_st_xsq_mosi : IN  t_mem_mosi;  
+    ram_st_xsq_miso : OUT t_mem_miso
+  );
+END st_xsq_arr;
+
+ARCHITECTURE str OF st_xsq_arr IS
+
+  CONSTANT c_xsq : NATURAL := g_nof_signal_inputs * g_nof_signal_inputs;
+  CONSTANT c_nof_statistics : NATURAL := g_nof_crosslets * c_xsq;
+  CONSTANT c_nof_word     : NATURAL := g_stat_data_sz*c_nof_statistics*c_nof_complex;
+  CONSTANT c_nof_word_w   : NATURAL := ceil_log2(c_nof_word);
+  
+  SIGNAL ram_st_xsq_mosi_arr : t_mem_mosi_arr(g_nof_streams-1 DOWNTO 0) := (OTHERS => c_mem_mosi_rst);
+  SIGNAL ram_st_xsq_miso_arr : t_mem_miso_arr(g_nof_streams-1 DOWNTO 0) := (OTHERS => c_mem_miso_rst);
+    
+BEGIN
+  -- st_xsq instances
+  gen_xsq : FOR I IN 0 TO g_nof_streams-1 GENERATE
+    st_xsq : ENTITY work.st_xsq 
+    GENERIC MAP (
+      g_nof_signal_inputs => g_nof_signal_inputs,
+      g_nof_crosslets     => g_nof_crosslets,
+      g_in_data_w         => g_in_data_w,   
+      g_stat_data_w       => g_stat_data_w,      
+      g_stat_data_sz      => g_stat_data_sz     
+    )
+    PORT MAP (
+  
+      mm_rst => mm_rst, 
+      mm_clk => mm_clk, 
+      dp_rst => dp_rst, 
+      dp_clk => dp_clk, 
+               
+      -- Streaming    
+      in_a => in_a_arr(I),  
+      in_b => in_b_arr(I),  
+      
+      -- Memory Mapped
+      ram_st_xsq_mosi => ram_st_xsq_mosi_arr(I), 
+      ram_st_xsq_miso => ram_st_xsq_miso_arr(I) 
+    );
+  END GENERATE;
+  
+  ---------------------------------------------------------------
+  -- COMBINE MEMORY MAPPED INTERFACES
+  ---------------------------------------------------------------
+  -- Combine the internal array of mm interfaces.
+  u_mem_mux_select : entity common_lib.common_mem_mux
+  generic map (    
+    g_nof_mosi    => g_nof_streams,
+    g_mult_addr_w => c_nof_word_w
+  )
+  port map (
+    mosi     => ram_st_xsq_mosi,
+    miso     => ram_st_xsq_miso,
+    mosi_arr => ram_st_xsq_mosi_arr,
+    miso_arr => ram_st_xsq_miso_arr
+  );
+ 
+  
+END str;
diff --git a/libraries/dsp/st/tb/vhdl/tb_st_pkg.vhd b/libraries/dsp/st/tb/vhdl/tb_st_pkg.vhd
new file mode 100644
index 0000000000000000000000000000000000000000..7d6c17284ffb148c7eafc8f037b1f71e5d585752
--- /dev/null
+++ b/libraries/dsp/st/tb/vhdl/tb_st_pkg.vhd
@@ -0,0 +1,54 @@
+-------------------------------------------------------------------------------
+--
+-- Copyright 2021
+-- ASTRON (Netherlands Institute for Radio Astronomy) <http://www.astron.nl/>
+-- P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
+--
+-- Licensed under the Apache License, Version 2.0 (the "License");
+-- you may not use this file except in compliance with the License.
+-- You may obtain a copy of the License at
+--
+--     http://www.apache.org/licenses/LICENSE-2.0
+--
+-- Unless required by applicable law or agreed to in writing, software
+-- distributed under the License is distributed on an "AS IS" BASIS,
+-- WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+-- See the License for the specific language governing permissions and
+-- limitations under the License.
+--
+-------------------------------------------------------------------------------
+-- Author: R. vd Walle
+-- Purpose:
+-- Functions used in st_lib 
+LIBRARY IEEE, common_lib;
+USE IEEE.std_logic_1164.ALL;
+USE IEEE.numeric_std.ALL;
+USE common_lib.common_pkg.ALL;
+
+PACKAGE tb_st_pkg IS 
+
+  FUNCTION func_st_calculate_expected_xsq(a_re, a_im, b_re, b_im : t_integer_arr; N_crosslets, N_int : NATURAL) RETURN t_integer_arr; 
+ 
+END tb_st_pkg;
+
+PACKAGE BODY tb_st_pkg IS 
+
+  FUNCTION func_st_calculate_expected_xsq(a_re, a_im, b_re, b_im : t_integer_arr; N_crosslets, N_int : NATURAL) RETURN t_integer_arr IS
+    CONSTANT c_N_s : NATURAL := a_re'LENGTH/N_crosslets;
+    CONSTANT c_xsq : NATURAL := c_N_s * c_N_s;
+    CONSTANT c_N_stat : NATURAL := N_crosslets * c_xsq;
+    VARIABLE v_exp_xsq : t_integer_arr(0 TO c_nof_complex * c_N_stat-1);
+  BEGIN
+
+    FOR N IN 0 TO N_crosslets-1 LOOP
+      FOR I IN 0 TO c_N_s-1 LOOP
+        FOR J IN 0 TO c_N_s-1 LOOP
+          v_exp_xsq(c_nof_complex * (N * c_xsq + I * c_N_s + J)    ) := N_int * COMPLEX_MULT_REAL(a_re(N * c_N_s + I), a_im(N * c_N_s + I), b_re(N * c_N_s + J), -1 * b_im(N * c_N_s + J));
+          v_exp_xsq(c_nof_complex * (N * c_xsq + I * c_N_s + J) + 1) := N_int * COMPLEX_MULT_IMAG(a_re(N * c_N_s + I), a_im(N * c_N_s + I), b_re(N * c_N_s + J), -1 * b_im(N * c_N_s + J));
+        END LOOP;
+      END LOOP;
+    END LOOP;
+    RETURN v_exp_xsq;
+  END;
+ 
+END tb_st_pkg;
diff --git a/libraries/dsp/st/tb/vhdl/tb_st_xsq.vhd b/libraries/dsp/st/tb/vhdl/tb_st_xsq.vhd
new file mode 100644
index 0000000000000000000000000000000000000000..8bb8ecc75340026968bfc873f4583172de7391a9
--- /dev/null
+++ b/libraries/dsp/st/tb/vhdl/tb_st_xsq.vhd
@@ -0,0 +1,232 @@
+-------------------------------------------------------------------------------
+--
+-- Copyright 2021
+-- ASTRON (Netherlands Institute for Radio Astronomy) <http://www.astron.nl/>
+-- P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
+--
+-- Licensed under the Apache License, Version 2.0 (the "License");
+-- you may not use this file except in compliance with the License.
+-- You may obtain a copy of the License at
+--
+--     http://www.apache.org/licenses/LICENSE-2.0
+--
+-- Unless required by applicable law or agreed to in writing, software
+-- distributed under the License is distributed on an "AS IS" BASIS,
+-- WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+-- See the License for the specific language governing permissions and
+-- limitations under the License.
+--
+-------------------------------------------------------------------------------
+
+-------------------------------------------------------------------------------
+-- Author : R vd Walle
+-- Purpose:  Testbench for the st_xsq unit. 
+--
+-- Usage in non-auto-mode (c_modelsim_start = 0 in python):
+--   > as 5
+--   > run -all
+-- Description: 
+-- The tb generates random data to feed into st_xsq. The output is compared to
+-- a pre-calculated expected array of xsq values.
+-- Remark:
+-- . More detail can be found in:
+--   https://support.astron.nl/confluence/display/L2M/L5+SDPFW+Design+Document%3A+Subband+Correlator
+
+LIBRARY IEEE, common_lib, mm_lib, diag_lib, dp_lib;
+USE IEEE.std_logic_1164.ALL;
+USE IEEE.numeric_std.ALL;
+USE common_lib.common_pkg.ALL;
+USE common_lib.common_math_pkg.ALL;
+USE common_lib.common_mem_pkg.ALL;
+USE common_lib.common_str_pkg.ALL;
+USE common_lib.tb_common_pkg.ALL;
+USE common_lib.tb_common_mem_pkg.ALL;
+USE common_lib.common_lfsr_sequences_pkg.ALL;
+USE mm_lib.mm_file_unb_pkg.ALL;
+USE mm_lib.mm_file_pkg.ALL;
+USE dp_lib.dp_stream_pkg.ALL; 
+USE diag_lib.diag_pkg.ALL; 
+USE dp_lib.tb_dp_pkg.ALL;
+USE work.tb_st_pkg.ALL;
+
+ENTITY tb_st_xsq IS 
+  GENERIC(
+    g_nof_crosslets      : NATURAL := 2; 
+    g_nof_signal_inputs  : NATURAL := 12;  
+    g_in_data_w          : NATURAL := 16;
+    g_nof_sync           : NATURAL := 3;
+    g_stat_data_w        : NATURAL := 64;  -- statistics accumulator width
+    g_stat_data_sz       : NATURAL := 2;   -- statistics word width >= statistics accumulator width and fit in a power of 2 multiple 32b MM words
+    g_nof_block_per_sync : NATURAL := 5;
+    g_nof_clk_per_blk    : NATURAL := 1024
+  );
+END tb_st_xsq;
+
+ARCHITECTURE tb OF tb_st_xsq IS
+  
+  CONSTANT c_sim            : BOOLEAN := TRUE;
+  CONSTANT c_rl             : NATURAL := 1;
+  CONSTANT c_block_size     : NATURAL := g_nof_crosslets * g_nof_signal_inputs;
+  CONSTANT c_random_data_w  : NATURAL := 8;
+  CONSTANT c_xsq            : NATURAL := g_nof_signal_inputs * g_nof_signal_inputs;
+  CONSTANT c_nof_statistics : NATURAL := g_nof_crosslets * c_xsq;
+  CONSTANT c_gap_size       : NATURAL := g_nof_clk_per_blk - c_nof_statistics;
+
+  CONSTANT c_random_in_a_re : t_integer_arr(0 TO c_block_size-1) := common_math_create_random_arr(c_block_size, c_random_data_w, 100);
+  CONSTANT c_random_in_a_im : t_integer_arr(0 TO c_block_size-1) := common_math_create_random_arr(c_block_size, c_random_data_w, 101);
+  CONSTANT c_random_in_b_re : t_integer_arr(0 TO c_block_size-1) := common_math_create_random_arr(c_block_size, c_random_data_w, 102);
+  CONSTANT c_random_in_b_im : t_integer_arr(0 TO c_block_size-1) := common_math_create_random_arr(c_block_size, c_random_data_w, 103);
+
+  CONSTANT c_expected_xsq   : t_integer_arr(0 TO c_nof_statistics * c_nof_complex-1) := func_st_calculate_expected_xsq(c_random_in_a_re, c_random_in_a_im, c_random_in_b_re, c_random_in_b_im, g_nof_crosslets, g_nof_block_per_sync);
+  ----------------------------------------------------------------------------
+  -- Clocks and resets
+  ----------------------------------------------------------------------------   
+  CONSTANT c_mm_clk_period  : TIME := 100 ps;
+  CONSTANT c_dp_clk_period  : TIME := 5 ns;
+  CONSTANT c_dp_pps_period  : NATURAL := 64;
+
+  SIGNAL tb_end             : STD_LOGIC;
+  SIGNAL dp_pps             : STD_LOGIC;
+
+  SIGNAL mm_rst             : STD_LOGIC := '1';
+  SIGNAL mm_clk             : STD_LOGIC := '1';
+
+  SIGNAL dp_rst             : STD_LOGIC;
+  SIGNAL dp_clk             : STD_LOGIC := '1';
+
+  SIGNAL st_en              : STD_LOGIC := '1';
+  SIGNAL st_siso            : t_dp_siso := c_dp_siso_rdy;
+  SIGNAL st_sosi            : t_dp_sosi := c_dp_sosi_rst;
+  SIGNAL in_sosi_a          : t_dp_sosi := c_dp_sosi_rst;
+  SIGNAL in_sosi_b          : t_dp_sosi := c_dp_sosi_rst;
+
+  ----------------------------------------------------------------------------
+  -- MM buses
+  ----------------------------------------------------------------------------                                         
+  SIGNAL ram_st_xsq_mosi    : t_mem_mosi := c_mem_mosi_rst;
+  SIGNAL ram_st_xsq_miso    : t_mem_miso := c_mem_miso_rst;
+  SIGNAL st_xsq_out_arr     : t_slv_32_arr(0 TO c_nof_statistics * c_nof_complex * g_stat_data_sz-1) := (OTHERS => (OTHERS => '0'));
+BEGIN
+
+  ----------------------------------------------------------------------------
+  -- Clock and reset generation
+  ----------------------------------------------------------------------------
+  mm_clk <= NOT mm_clk OR tb_end AFTER c_mm_clk_period/2;
+  mm_rst <= '1', '0' AFTER c_mm_clk_period*5;
+
+  dp_clk <= NOT dp_clk OR tb_end AFTER c_dp_clk_period/2;
+  dp_rst <= '1', '0' AFTER c_dp_clk_period*5;
+
+  ------------------------------------------------------------------------------
+  -- MM Stimuli
+  ------------------------------------------------------------------------------
+  p_mm_stimuli : PROCESS
+  BEGIN
+    FOR I IN 0 TO g_nof_sync -1 LOOP
+      proc_common_wait_until_lo_hi(dp_clk, st_sosi.sync);
+    END LOOP;
+    proc_common_wait_some_cycles(dp_clk, g_nof_clk_per_blk);
+    proc_common_wait_some_cycles(dp_clk, 20);
+    -- read ram
+    proc_mem_read_ram(0, c_nof_statistics * c_nof_complex * g_stat_data_sz, mm_clk, ram_st_xsq_mosi, ram_st_xsq_miso, st_xsq_out_arr);  
+    WAIT;
+  END PROCESS;
+
+  ------------------------------------------------------------------------------
+  -- Data blocks
+  ------------------------------------------------------------------------------
+  p_st_stimuli : PROCESS
+    VARIABLE v_re  : NATURAL := 0;
+    VARIABLE v_im  : NATURAL := 0;
+  BEGIN  
+    tb_end <= '0';
+    st_sosi <= c_dp_sosi_rst;
+    proc_common_wait_until_low(dp_clk, dp_rst);
+  
+    -- Run some sync intervals with DSP counter data for the real and imag fields
+    WAIT UNTIL rising_edge(dp_clk);
+    proc_common_wait_some_cycles(dp_clk, 7);
+    FOR I IN 0 TO g_nof_sync-1 LOOP
+      proc_dp_gen_block_data(c_rl, FALSE, g_in_data_w, g_in_data_w, 0, v_re, v_im, c_nof_statistics, 0, 0, '1', "0", dp_clk, st_en, st_siso, st_sosi);     -- next sync
+      st_sosi <= c_dp_sosi_rst;
+      proc_common_wait_some_cycles(dp_clk, c_gap_size);
+      FOR J IN 0 TO g_nof_block_per_sync-2 LOOP  -- provide sop and eop for block reference
+        proc_dp_gen_block_data(c_rl, FALSE, g_in_data_w, g_in_data_w, 0, v_re, v_im, c_nof_statistics, 0, 0, '0', "0", dp_clk, st_en, st_siso, st_sosi);   -- no sync
+        st_sosi <= c_dp_sosi_rst;
+        proc_common_wait_some_cycles(dp_clk, c_gap_size);
+      END LOOP;
+    END LOOP;
+    st_sosi <= c_dp_sosi_rst;
+    proc_common_wait_some_cycles(dp_clk, 100);
+    tb_end <= '1';
+    WAIT;
+  END PROCESS;
+
+  in_sosi_a.sync   <= st_sosi.sync;
+  in_sosi_b.sync   <= st_sosi.sync;
+  in_sosi_a.sop    <= st_sosi.sop;
+  in_sosi_b.sop    <= st_sosi.sop;
+  in_sosi_a.eop    <= st_sosi.eop;
+  in_sosi_b.eop    <= st_sosi.eop;
+  in_sosi_a.valid  <= st_sosi.valid;
+  in_sosi_b.valid  <= st_sosi.valid;
+
+  p_random_stimuli : PROCESS
+  BEGIN
+    FOR S IN 0 TO g_nof_sync * g_nof_block_per_sync -1 LOOP
+      WAIT UNTIL rising_edge(st_sosi.sop);
+      FOR N IN 0 TO g_nof_crosslets -1 LOOP
+        FOR I IN 0 TO g_nof_signal_inputs -1 LOOP
+          FOR J IN 0 TO g_nof_signal_inputs -1 LOOP      
+            in_sosi_a.re <= TO_DP_DSP_DATA(c_random_in_a_re(N * g_nof_signal_inputs + I));
+            in_sosi_a.im <= TO_DP_DSP_DATA(c_random_in_a_im(N * g_nof_signal_inputs + I));
+            in_sosi_b.re <= TO_DP_DSP_DATA(c_random_in_b_re(N * g_nof_signal_inputs + J));
+            in_sosi_b.im <= TO_DP_DSP_DATA(c_random_in_b_im(N * g_nof_signal_inputs + J));
+            proc_common_wait_some_cycles(dp_clk, 1);
+          END LOOP;
+        END LOOP;
+      END LOOP;
+    END LOOP;
+    WAIT;
+  END PROCESS;
+
+  ------------------------------------------------------------------------------
+  -- Verification 
+  ------------------------------------------------------------------------------
+  p_verify : PROCESS
+  BEGIN
+    proc_common_wait_until_high(mm_clk, ram_st_xsq_miso.rdval);
+    proc_common_wait_some_cycles(mm_clk, 2* c_nof_statistics * c_nof_complex * g_stat_data_sz + 10);
+    FOR I IN 0 TO c_nof_statistics * c_nof_complex -1 LOOP
+      ASSERT TO_SINT(st_xsq_out_arr(I*2)) = c_expected_xsq(I) REPORT "WRONG XSQ DATA" SEVERITY ERROR; -- Only read low part of statistic
+    END LOOP;
+    WAIT;
+  END PROCESS;
+
+  ----------------------------------------------------------------------------
+  -- DUT: Device Under Test
+  ---------------------------------------------------------------------------- 
+  u_dut : ENTITY work.st_xsq
+  GENERIC MAP(
+    g_nof_signal_inputs  => g_nof_signal_inputs,      
+    g_nof_crosslets      => g_nof_crosslets,    
+    g_in_data_w          => g_in_data_w,     
+    g_stat_data_w        => g_stat_data_w,   
+    g_stat_data_sz       => g_stat_data_sz  
+  )
+  PORT MAP(
+    mm_rst           =>  mm_rst,
+    mm_clk           =>  mm_clk,
+    dp_rst           =>  dp_rst,
+    dp_clk           =>  dp_clk,
+    
+    -- Streaming
+    in_a             => in_sosi_a, 
+    in_b             => in_sosi_b, 
+    
+    -- Memory Mapped
+    ram_st_xsq_mosi  => ram_st_xsq_mosi,
+    ram_st_xsq_miso  => ram_st_xsq_miso
+  );
+
+END tb;
diff --git a/libraries/dsp/st/tb/vhdl/tb_tb_st_xsq.vhd b/libraries/dsp/st/tb/vhdl/tb_tb_st_xsq.vhd
new file mode 100644
index 0000000000000000000000000000000000000000..ab1b84dfa8f35bea402a0d2be7c582845ebb8731
--- /dev/null
+++ b/libraries/dsp/st/tb/vhdl/tb_tb_st_xsq.vhd
@@ -0,0 +1,56 @@
+-------------------------------------------------------------------------------
+--
+-- Copyright 2021
+-- ASTRON (Netherlands Institute for Radio Astronomy) <http://www.astron.nl/>
+-- P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
+--
+-- Licensed under the Apache License, Version 2.0 (the "License");
+-- you may not use this file except in compliance with the License.
+-- You may obtain a copy of the License at
+--
+--     http://www.apache.org/licenses/LICENSE-2.0
+--
+-- Unless required by applicable law or agreed to in writing, software
+-- distributed under the License is distributed on an "AS IS" BASIS,
+-- WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+-- See the License for the specific language governing permissions and
+-- limitations under the License.
+--
+-------------------------------------------------------------------------------
+
+-------------------------------------------------------------------------------
+-- Author : R vd Walle
+-- Purpose: Test multiple instances of tb_st_xsq
+-- Usage:
+-- > as 10
+-- > run -all
+--
+-- Description: See tb_st_xsq
+
+LIBRARY IEEE;
+USE IEEE.std_logic_1164.ALL;
+
+ENTITY tb_tb_st_xsq IS
+END tb_tb_st_xsq;
+
+ARCHITECTURE tb OF tb_tb_st_xsq IS
+
+  CONSTANT c_nof_sync       : NATURAL := 3;
+  CONSTANT c_dsp_data_w     : NATURAL := 16;
+  SIGNAL tb_end : STD_LOGIC := '0';  -- declare tb_end to avoid 'No objects found' error on 'when -label tb_end'
+  
+BEGIN
+--  GENERICS:
+--    g_nof_crosslets      : NATURAL := 2; 
+--    g_nof_signal_inputs  : NATURAL := 12;  
+--    g_in_data_w          : NATURAL := 16;
+--    g_nof_sync           : NATURAL := 3;
+--    g_stat_data_w        : NATURAL := 64;  -- statistics accumulator width
+--    g_stat_data_sz       : NATURAL := 2;   -- statistics word width >= statistics accumulator width and fit in a power of 2 multiple 32b MM words
+--    g_nof_block_per_sync : NATURAL := 5;
+--    g_nof_clk_per_blk    : NATURAL := 1024
+
+  u_sdp : ENTITY work.tb_st_xsq GENERIC MAP (1, 12, c_dsp_data_w, c_nof_sync, 64, 2, 5, 1024);
+  u_max : ENTITY work.tb_st_xsq GENERIC MAP (16, 8, c_dsp_data_w, c_nof_sync, 64, 2, 5, 1024); -- g_nof_crosslets * g_nof_signal_inputs**2 = 16 * 8 * 8 = 1024 = g_nof_clk_per_blk 
+  
+END tb;