Subversion Repositories pulse_processing_algorithm

-----------------------------------------------------------------------------------------------
--
--    Copyright (C) 2011 Peter Lemmens, PANDA collaboration
--		p.j.j.lemmens@rug.nl
--    http://www-panda.gsi.de
--
--    As a reference, please use:
--    E. Guliyev, M. Kavatsyuk, P.J.J. Lemmens, G. Tambave, H. Loehner,
--    "VHDL Implementation of Feature-Extraction Algorithm for the PANDA Electromagnetic Calorimeter"
--    Nuclear Inst. and Methods in Physics Research, A ....
--
--
--    This program is free software; you can redistribute it and/or modify
--    it under the terms of the GNU Lesser General Public License as published by
--    the Free Software Foundation; either version 3 of the License, or
--    (at your option) any later version.
--
--    This program is distributed in the hope that it will be useful,
--    but WITHOUT ANY WARRANTY; without even the implied warranty of
--    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
--    GNU Lesser General Public License for more details.
--
--    You should have received a copy of the GNU General Public License
--    along with this program; if not, write to the Free Software
--    Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111 USA
--
-----------------------------------------------------------------------------------------------
-----------------------------------------------------------------------------------------------
-- Company:			KVI (Kernfysisch Versneller Instituut  -- Groningen, The Netherlands	
-- Author:			P.J.J. Lemmens
-- Design Name:	mwd_programmable
-- Module Name:	MWD_programmable.vhd
-- Description:	moving-window deconvolution algorithm (MWD) with programmable 
--						filter-length (power of 2) and decay-correction parameter
--						Intended for use with 'slow' logarithmic-decaying pulses that 
--						occur for instance when a particle/photon strikes a photo-diode
-----------------------------------------------------------------------------------------------
library IEEE;
use IEEE.STD_LOGIC_1164.ALL;
use IEEE.STD_LOGIC_ARITH.ALL;
use IEEE.STD_LOGIC_SIGNED.ALL;
 
entity mwd_programmable is
	generic(	MAX_MWD_PWR		:	natural);
   Port (rst					: in  STD_LOGIC;
			clk					: in  STD_LOGIC;
			enable				: in  STD_LOGIC;
			program				: in  std_logic;
			mwd_pwr_in			: in  std_logic_vector(7 downto 0);
			data_in				: in  STD_LOGIC_VECTOR;
			correction_in		: in  STD_LOGIC_VECTOR;
			data_out				: out STD_LOGIC_VECTOR
		);
end mwd_programmable;
 
architecture Behavioral of mwd_programmable is
 
	constant	WIDTH			: natural := data_in'length;
	constant	M_WIDTH		: natural := 18;							-- Multiplier width ; xilinx MULT18x18_SIO
 
	component window_subtractor_programmable
		generic(	FORWARD		:	boolean;
					MAX_MWD_PWR	:	natural);
		Port (rst				: in  STD_LOGIC;
				clk				: in  STD_LOGIC;
				enable			: IN	STD_LOGIC ;
				program			: in  STD_LOGIC;
				window_pwr_in	: in  STD_LOGIC_VECTOR(7 downto 0);
				data_in			: in  STD_LOGIC_VECTOR;
				data_out			: out STD_LOGIC_VECTOR
				);
	end component;
 
	component moving_average_programmable
		generic(MEM_PWR		: natural);
		port (rst				: in	std_logic;
				clk				: in	std_logic;
				enable			: in  std_logic;
				program			: in  std_logic;
				avg_pwr_in		: in  std_logic_vector(7 downto 0);
				data_in			: in  std_logic_vector;
				data_out			: out std_logic_vector
			);
	end component;
 
----------------------------------------------------------------------------------------------------------
--The calculations of and with correction factor introduce interesting problems combined with
--the fact that the dedicated hardware multipliers in Xilinx Spartan III are 18/36bit wide(fixed)
--
--(1) As we have to multiply with a sum over (2^WINDOW_PWR) samples, the sum should be 
--(WIDTH+WINDOW_PWR)-bits wide. However if this width exceeds 18bit we will have to truncate the result,
--effectively dividing by (2^(lost#bits)). We can identify 3 situations where (WIDTH+WINDOW_PWR)<18,
--(WIDTH+WINDOW_PWR)=18 and (WIDTH+WINDOW_PWR)>18 that require different approaches making a generic and
--scalable solution all but impossible.
--
--(2) The correction factor on the other hand is a real number <1 which we want to convert to a usefull
--integer value in order to be able to do a cost-effective multiplication. In the given example/situation
--correction_factor = 0.69314718/tau = 0.69314718 / 1250 ~ 0.554517e-3. This allows for a multiplication
--by 2^28 while still remaining within 18bits (as 262143 is the max int in 18bit); result: 148852
--
----------------------------------------------------------------------------------------------------------
--
--Example for 16-bit data and a WINDOW_PWR of 6 (64 sample window):
--** ((WIDTH+WINDOW_PWR)= 16 + 6 = 22bit, fitting into 18bit means a loss of 4-bit equals a division by 16 !!
--** Correction factor multiplication (2) gives us a gain of 28-bit, thus effectively the PRODUCT will be 2^24 too high
--
--************************************************************
--**** the lower 24 bits should therefore be discarded !! ****
--************************************************************
--
----------------------------------------------------------------------------------------------------------
 
 
--	constant Fsmp	: real := 50.0e6;										-- Sample frequency 50MHz
--	constant tau	: real := (0.0000125	* Fsmp);							-- calculate time-constant
--	constant ln2	: real := 0.69314718;									-- natural log of 2
--	constant correction_factor	: integer range 0 to (2**18 - 1) := integer((ln2/tau) * real(2**24));	-- decay correction weight*(2^24)
--	constant correction_factor	: std_logic_vector(17  downto 0) := conv_std_logic_vector(signed(148852), 18);
--	constant CORRECTION_FACTOR	: std_logic_vector(17  downto 0) := conv_std_logic_vector(74426, 18); -- = (ln(2)/tau) * 2^26)
	-- this should give a integer value of correction_factor=18607 which still fits
	-- within the 18bit (2^18=262144)... 2^27 is max !!
 
	signal rst_S				: std_logic := '1';
	signal clk_S				: std_logic := '0';
	signal enable_S			: std_logic := '0';
	signal data_in_S			: std_logic_vector(WIDTH - 1 downto 0) := (others => '0');
	signal program_S			: std_logic := '0';
	signal mwd_pwr_S			: std_logic_vector(7 downto 0) 				:= x"05";
	signal diff_data_S		: std_logic_vector(WIDTH - 1 downto 0)		:= (others => '0');
	signal decay_corr_S		: std_logic_vector(WIDTH - 1 downto 0)	:= (others => '0');
	signal average_S			: std_logic_vector(WIDTH - 1 downto 0)		:= (others => '0');
	signal mult_sum_S			: std_logic_vector(M_WIDTH - 1 downto 0)	:= (others => '0');
	signal factor_S			: std_logic_vector(M_WIDTH - 1 downto 0)	:= conv_std_logic_vector(37213, M_WIDTH);
	signal product36			: std_logic_vector(35 downto 0)				:= (others => '0');
	signal data_out_S			: std_logic_vector(WIDTH - 1 downto 0)		:= (others => '0');
 
-----------------------------------------------------------------------
 
	begin
 
		rdif : window_subtractor_programmable
			generic map(FORWARD			=> true,
							MAX_MWD_PWR		=> MAX_MWD_PWR)
			Port map(	rst				=> rst_S,
							clk				=> clk_S,
							enable 			=>	enable_S,
							program			=> program_S,
							window_pwr_in	=> mwd_pwr_S,
							data_in			=> data_in_S,
							data_out 		=> diff_data_S	
						);
 
		mwd_sample_sum : moving_average_programmable
			generic map(MEM_PWR	=>	MAX_MWD_PWR)		
			port map(	rst			=> rst_S,
							clk			=>	clk_S,
							enable		=>	enable_S,
							program		=> program_S,
							avg_pwr_in	=> mwd_pwr_S,
							data_in		=>	data_in_S,
							data_out		=> average_S
						);
 
 
		rst_S				<= rst;
		clk_S				<= clk;
		enable_S			<= enable;
					mwd_pwr_S	<= mwd_pwr_in;
					factor_S		<=	conv_std_logic_vector(conv_integer(correction_in), factor_S'length);
		program_S		<= program;		-- allready synced and single clock by flowcontrol
		data_in_S		<= data_in;
		data_out			<= data_out_S;
 
--	program_proc : process (clk_S)
--		begin
--		if rising_edge(clk_S) then
----			program_S		<= program;
--			if (rst_S = '1') then
--				mwd_pwr_S	<= x"05";
--			else
--				if (program = '1') then
--					if (mwd_pwr_in >= 2) and (mwd_pwr_in <= 6) then
--						mwd_pwr_S	<= mwd_pwr_in;
--					else
--						mwd_pwr_S	<= x"05";
--					end if;
--					factor_S		<=	conv_std_logic_vector(conv_integer(correction_in), factor_S'length);
--				end if;
--			end if;
--		end if;
--	end process;
 
-- all averages should be multiplied by  2^mwd_pwr (we really need a sum) but this generates overflows
-- instead we take the average and map it on 18 bit then we choose the fitting correction factor
 
		mult_sum_S	<= conv_std_logic_vector(conv_integer(average_S), M_WIDTH);
 
		--	READ explanation above  !!
		decay_corr_S(15 downto 0) 						<= product36(35 downto 20);		-- WATCH IT !! 36bits !! not 32 !! scale it down 24 bits
		decay_corr_S(decay_corr_S'high downto 16) <= (others => '0');					--  
 
		decay_correction_proc : process(clk_S, rst_S)
		begin
			if (clk_S'event and clk_S = '1') then
				if (rst_S = '1') then
					product36		<= (others => '0');
					data_out_S		<= (others => '0');
				else
					if (enable_S = '1') then
						product36	<= mult_sum_S * factor_S;
						data_out_S	<= diff_data_S + decay_corr_S;
					end if;
				end if;
			end if;
		end process;
 
end Behavioral;