From 094f35f19ddf83aed14a093db6c1b187d2d65012 Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Wed, 19 Mar 2025 18:11:01 -0600
Subject: [PATCH 01/20] add stim demultiplex rhs intan

---
 neo/rawio/intanrawio.py | 65 ++++++++++++++++++++++++++++++++++++++---
 1 file changed, 61 insertions(+), 4 deletions(-)

diff --git a/neo/rawio/intanrawio.py b/neo/rawio/intanrawio.py
index be2db6a59..222312cbd 100644
--- a/neo/rawio/intanrawio.py
+++ b/neo/rawio/intanrawio.py
@@ -443,6 +443,7 @@ def _get_analogsignal_chunk_header_attached(self, i_start, i_stop, stream_index,
 
         stream_name = self.header["signal_streams"][stream_index]["name"][:]
         stream_is_digital = stream_name in digital_stream_names
+        stream_is_stim = stream_name == "Stim channel"
 
         field_name = stream_name if stream_is_digital else channel_ids[0]
 
@@ -462,7 +463,18 @@ def _get_analogsignal_chunk_header_attached(self, i_start, i_stop, stream_index,
             sl1 = sl0 + (i_stop - i_start)
 
         # For all streams raw_data is a structured memmap with a field for each channel_id
-        if not stream_is_digital:
+        if stream_is_stim:
+            # For stim data, we need to extract the raw data first, then demultiplex it
+            stim_data = np.zeros((i_stop - i_start, len(channel_ids)), dtype=dtype)
+            for chunk_index, channel_id in enumerate(channel_ids):
+                data_chan = self._raw_data[channel_id]
+                if multiple_samples_per_block:
+                    stim_data[:, chunk_index] = data_chan[block_start:block_stop].flatten()[sl0:sl1]
+                else:
+                    stim_data[:, chunk_index] = data_chan[i_start:i_stop]
+            # Now demultiplex the stim data
+            sigs_chunk = self._demultiplex_stim_data(stim_data, 0, stim_data.shape[0])
+        elif not stream_is_digital:
             sigs_chunk = np.zeros((i_stop - i_start, len(channel_ids)), dtype=dtype)
 
             for chunk_index, channel_id in enumerate(channel_ids):
@@ -480,6 +492,8 @@ def _get_analogsignal_chunk_one_file_per_channel(self, i_start, i_stop, stream_i
 
         stream_name = self.header["signal_streams"][stream_index]["name"][:]
         signal_data_memmap_list = self._raw_data[stream_name]
+        stream_is_stim = stream_name == "Stim channel"
+        
         channel_indexes_are_slice = isinstance(channel_indexes, slice)
         if channel_indexes_are_slice:
             num_channels = len(signal_data_memmap_list)
@@ -496,6 +510,10 @@ def _get_analogsignal_chunk_one_file_per_channel(self, i_start, i_stop, stream_i
         for chunk_index, channel_index in enumerate(channel_indexes):
             channel_memmap = signal_data_memmap_list[channel_index]
             sigs_chunk[:, chunk_index] = channel_memmap[i_start:i_stop]
+            
+        # If this is stim data, we need to demultiplex it
+        if stream_is_stim:
+            sigs_chunk = self._demultiplex_stim_data(sigs_chunk, 0, sigs_chunk.shape[0])
 
         return sigs_chunk
 
@@ -505,6 +523,8 @@ def _get_analogsignal_chunk_one_file_per_signal(self, i_start, i_stop, stream_in
         raw_data = self._raw_data[stream_name]
 
         stream_is_digital = stream_name in digital_stream_names
+        stream_is_stim = stream_name == "Stim channel"
+        
         if stream_is_digital:
             stream_id = self.header["signal_streams"][stream_index]["id"]
             mask = self.header["signal_channels"]["stream_id"] == stream_id
@@ -512,7 +532,8 @@ def _get_analogsignal_chunk_one_file_per_signal(self, i_start, i_stop, stream_in
             channel_ids = signal_channels["id"][channel_indexes]
 
             output = self._demultiplex_digital_data(raw_data, channel_ids, i_start, i_stop)
-
+        elif stream_is_stim:
+            output = self._demultiplex_stim_data(raw_data, i_start, i_stop)
         else:
             output = raw_data[i_start:i_stop, channel_indexes]
 
@@ -530,6 +551,42 @@ def _demultiplex_digital_data(self, raw_digital_data, channel_ids, i_start, i_st
             output[:, channel_index] = demultiplex_data[i_start:i_stop].flatten()
 
         return output
+        
+    def _demultiplex_stim_data(self, raw_stim_data, i_start, i_stop):
+        """
+        Demultiplexes the stim data stream.
+        
+        Parameters
+        ----------
+        raw_stim_data : ndarray
+            The raw stim data
+        i_start : int
+            Start index
+        i_stop : int
+            Stop index
+            
+        Returns
+        -------
+        output : ndarray
+            Demultiplexed stim data containing only the current values, preserving channel dimensions
+        """
+        # Get the relevant portion of the data
+        data = raw_stim_data[i_start:i_stop]
+        
+        # Extract current value (bits 0-8)
+        magnitude = np.bitwise_and(data, 0xFF)  # Extract lowest 8 bits
+        sign_bit = np.bitwise_and(np.right_shift(data, 8), 0x01)  # Extract 9th bit for sign
+        
+        # Apply sign to current values
+        current = np.where(sign_bit == 1, -magnitude, magnitude)
+        
+        # Note: If needed, other flag bits could be extracted as follows:
+        # compliance_flag = np.bitwise_and(np.right_shift(data, 15), 0x01).astype(bool)  # Bit 16 (MSB)
+        # charge_recovery_flag = np.bitwise_and(np.right_shift(data, 14), 0x01).astype(bool)  # Bit 15
+        # amp_settle_flag = np.bitwise_and(np.right_shift(data, 13), 0x01).astype(bool)  # Bit 14
+        # These could be returned as a structured array or dictionary if needed
+        
+        return current
 
     def get_intan_timestamps(self, i_start=None, i_stop=None):
         """
@@ -857,8 +914,8 @@ def read_rhs(filename, file_format: str):
             chan_info_stim["sampling_rate"] = sr
             # stim channel are complicated because they are coded
             # with bits, they do not fit the gain/offset rawio strategy
-            chan_info_stim["units"] = ""
-            chan_info_stim["gain"] = 1.0
+            chan_info_stim["units"] = "A"  # Amps
+            chan_info_stim["gain"] = global_info["stim_step_size"]
             chan_info_stim["offset"] = 0.0
             chan_info_stim["signal_type"] = 11  # put it in another group
             chan_info_stim["dtype"] = "uint16"

From 6dadc35bf58b2ecd6d6e09d985e598fe66323420 Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Wed, 19 Mar 2025 19:41:12 -0600
Subject: [PATCH 02/20] fix tests

---
 neo/rawio/intanrawio.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/neo/rawio/intanrawio.py b/neo/rawio/intanrawio.py
index 222312cbd..cd6c08519 100644
--- a/neo/rawio/intanrawio.py
+++ b/neo/rawio/intanrawio.py
@@ -534,6 +534,7 @@ def _get_analogsignal_chunk_one_file_per_signal(self, i_start, i_stop, stream_in
             output = self._demultiplex_digital_data(raw_data, channel_ids, i_start, i_stop)
         elif stream_is_stim:
             output = self._demultiplex_stim_data(raw_data, i_start, i_stop)
+            output = output[:, channel_indexes]
         else:
             output = raw_data[i_start:i_stop, channel_indexes]
 

From f1cfeb4884b7b94528ca27f17c1c5e8c8d3b4455 Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Wed, 19 Mar 2025 22:20:13 -0600
Subject: [PATCH 03/20] docstring and demultiplex to decode name change

---
 neo/rawio/intanrawio.py | 77 ++++++++++++++++++++++++++++++++++-------
 1 file changed, 64 insertions(+), 13 deletions(-)

diff --git a/neo/rawio/intanrawio.py b/neo/rawio/intanrawio.py
index cd6c08519..0b603fa76 100644
--- a/neo/rawio/intanrawio.py
+++ b/neo/rawio/intanrawio.py
@@ -473,7 +473,7 @@ def _get_analogsignal_chunk_header_attached(self, i_start, i_stop, stream_index,
                 else:
                     stim_data[:, chunk_index] = data_chan[i_start:i_stop]
             # Now demultiplex the stim data
-            sigs_chunk = self._demultiplex_stim_data(stim_data, 0, stim_data.shape[0])
+            sigs_chunk = self._decode_current_from_stim_data(stim_data, 0, stim_data.shape[0])
         elif not stream_is_digital:
             sigs_chunk = np.zeros((i_stop - i_start, len(channel_ids)), dtype=dtype)
 
@@ -513,7 +513,7 @@ def _get_analogsignal_chunk_one_file_per_channel(self, i_start, i_stop, stream_i
             
         # If this is stim data, we need to demultiplex it
         if stream_is_stim:
-            sigs_chunk = self._demultiplex_stim_data(sigs_chunk, 0, sigs_chunk.shape[0])
+            sigs_chunk = self._decode_current_from_stim_data(sigs_chunk, 0, sigs_chunk.shape[0])
 
         return sigs_chunk
 
@@ -533,7 +533,7 @@ def _get_analogsignal_chunk_one_file_per_signal(self, i_start, i_stop, stream_in
 
             output = self._demultiplex_digital_data(raw_data, channel_ids, i_start, i_stop)
         elif stream_is_stim:
-            output = self._demultiplex_stim_data(raw_data, i_start, i_stop)
+            output = self._decode_current_from_stim_data(raw_data, i_start, i_stop)
             output = output[:, channel_indexes]
         else:
             output = raw_data[i_start:i_stop, channel_indexes]
@@ -541,7 +541,39 @@ def _get_analogsignal_chunk_one_file_per_signal(self, i_start, i_stop, stream_in
         return output
 
     def _demultiplex_digital_data(self, raw_digital_data, channel_ids, i_start, i_stop):
-
+        """
+        Demultiplex digital data by extracting individual channel values from packed 16-bit format.
+        
+        According to the Intan format, digital input/output data is stored with all 16 channels
+        encoded bit-by-bit in each 16-bit word. This method extracts the specified digital channels
+        from the packed format into separate boolean arrays.
+        
+        Parameters
+        ----------
+        raw_digital_data : ndarray
+            Raw digital data in packed 16-bit format where each bit represents a different channel.
+        channel_ids : list or array
+            List of channel identifiers to extract. Each channel_id must correspond to a digital
+            input or output channel.
+        i_start : int
+            Starting sample index for demultiplexing.
+        i_stop : int
+            Ending sample index for demultiplexing (exclusive).
+            
+        Returns
+        -------
+        ndarray
+            Demultiplexed digital data with shape (i_stop-i_start, len(channel_ids)), 
+            containing boolean values for each requested channel.
+        
+        Notes
+        -----
+        In the Intan format, digital channels are packed into 16-bit words where each bit position
+        corresponds to a specific channel number. For example, with digital inputs 0, 4, and 5
+        set high and the rest low, the 16-bit word would be 2^0 + 2^4 + 2^5 = 1 + 16 + 32 = 49.
+        
+        The native_order property for each channel corresponds to its bit position in the packed word.
+        """
         dtype = np.uint16  # We fix this to match the memmap dtype
         output = np.zeros((i_stop - i_start, len(channel_ids)), dtype=dtype)
 
@@ -553,30 +585,49 @@ def _demultiplex_digital_data(self, raw_digital_data, channel_ids, i_start, i_st
 
         return output
         
-    def _demultiplex_stim_data(self, raw_stim_data, i_start, i_stop):
+    def _decode_current_from_stim_data(self, raw_stim_data, i_start, i_stop):
         """
-        Demultiplexes the stim data stream.
+        Demultiplex stimulation data by extracting current values from packed 16-bit format.
+        
+        According to the Intan RHS data format, stimulation current is stored in the lower 9 bits
+        of each 16-bit word: 8 bits for magnitude and 1 bit for sign. The upper bits contain
+        flags for compliance limit, charge recovery, and amplifier settle.
         
         Parameters
         ----------
         raw_stim_data : ndarray
-            The raw stim data
+            Raw stimulation data in packed 16-bit format.
         i_start : int
-            Start index
+            Starting sample index for demultiplexing.
         i_stop : int
-            Stop index
+            Ending sample index for demultiplexing (exclusive).
             
         Returns
         -------
-        output : ndarray
-            Demultiplexed stim data containing only the current values, preserving channel dimensions
+        ndarray
+            Demultiplexed stimulation current values in amperes, preserving the original
+            array dimensions. The output values need to be multiplied by the stim_step_size
+            parameter (from header) to obtain the actual current in amperes.
+        
+        Notes
+        -----
+        Bit structure of each 16-bit stimulation word:
+        - Bits 0-7: Current magnitude
+        - Bit 8: Sign bit (1 = negative current)
+        - Bits 9-13: Unused (always zero)
+        - Bit 14: Amplifier settle flag (1 = activated)
+        - Bit 15: Charge recovery flag (1 = activated) 
+        - Bit 16 (MSB): Compliance limit flag (1 = limit reached)
+        
+        The actual current value in amperes is obtained by multiplying the
+        output by the 'stim_step_size' parameter from the file header.
         """
         # Get the relevant portion of the data
         data = raw_stim_data[i_start:i_stop]
         
         # Extract current value (bits 0-8)
-        magnitude = np.bitwise_and(data, 0xFF)  # Extract lowest 8 bits
-        sign_bit = np.bitwise_and(np.right_shift(data, 8), 0x01)  # Extract 9th bit for sign
+        magnitude = np.bitwise_and(data, 0xFF)
+        sign_bit = np.bitwise_and(np.right_shift(data, 8), 0x01)  # Shift right by 8 bits to get the sign bit
         
         # Apply sign to current values
         current = np.where(sign_bit == 1, -magnitude, magnitude)

From be4a9e1630fd1fc16b270428133c776e80700d4f Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Thu, 20 Mar 2025 08:50:21 -0600
Subject: [PATCH 04/20] Update neo/rawio/intanrawio.py

Co-authored-by: Zach McKenzie <92116279+zm711@users.noreply.github.com>
---
 neo/rawio/intanrawio.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/neo/rawio/intanrawio.py b/neo/rawio/intanrawio.py
index 0b603fa76..a0ae65768 100644
--- a/neo/rawio/intanrawio.py
+++ b/neo/rawio/intanrawio.py
@@ -546,7 +546,7 @@ def _demultiplex_digital_data(self, raw_digital_data, channel_ids, i_start, i_st
         
         According to the Intan format, digital input/output data is stored with all 16 channels
         encoded bit-by-bit in each 16-bit word. This method extracts the specified digital channels
-        from the packed format into separate boolean arrays.
+        from the packed format into separate uint16 arrays of 0 and 1.
         
         Parameters
         ----------

From 2900bd34fed796aeb5ff55c4c1c9e366da283759 Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Thu, 20 Mar 2025 08:52:23 -0600
Subject: [PATCH 05/20] Update neo/rawio/intanrawio.py

Co-authored-by: Zach McKenzie <92116279+zm711@users.noreply.github.com>
---
 neo/rawio/intanrawio.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/neo/rawio/intanrawio.py b/neo/rawio/intanrawio.py
index a0ae65768..3a7c02247 100644
--- a/neo/rawio/intanrawio.py
+++ b/neo/rawio/intanrawio.py
@@ -564,7 +564,7 @@ def _demultiplex_digital_data(self, raw_digital_data, channel_ids, i_start, i_st
         -------
         ndarray
             Demultiplexed digital data with shape (i_stop-i_start, len(channel_ids)), 
-            containing boolean values for each requested channel.
+            containing 0 or 1 values for each requested channel.
         
         Notes
         -----

From f34680d3ec75685417eee88e086facf2582fae07 Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Thu, 20 Mar 2025 08:59:58 -0600
Subject: [PATCH 06/20] Update neo/rawio/intanrawio.py

Co-authored-by: Zach McKenzie <92116279+zm711@users.noreply.github.com>
---
 neo/rawio/intanrawio.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/neo/rawio/intanrawio.py b/neo/rawio/intanrawio.py
index 3a7c02247..b8286bf28 100644
--- a/neo/rawio/intanrawio.py
+++ b/neo/rawio/intanrawio.py
@@ -620,7 +620,8 @@ def _decode_current_from_stim_data(self, raw_stim_data, i_start, i_stop):
         - Bit 16 (MSB): Compliance limit flag (1 = limit reached)
         
         The actual current value in amperes is obtained by multiplying the
-        output by the 'stim_step_size' parameter from the file header.
+        output by the 'stim_step_size' parameter from the file header. These scaled values can be 
+        obtained with the `rescale_signal_raw_to_float` function.
         """
         # Get the relevant portion of the data
         data = raw_stim_data[i_start:i_stop]

From ea95545a18f1f81b12a7b6194de95d2e8bbbb450 Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Fri, 21 Mar 2025 11:51:01 -0600
Subject: [PATCH 07/20] fix

---
 neo/rawio/intanrawio.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/neo/rawio/intanrawio.py b/neo/rawio/intanrawio.py
index b8286bf28..4defea300 100644
--- a/neo/rawio/intanrawio.py
+++ b/neo/rawio/intanrawio.py
@@ -631,6 +631,9 @@ def _decode_current_from_stim_data(self, raw_stim_data, i_start, i_stop):
         sign_bit = np.bitwise_and(np.right_shift(data, 8), 0x01)  # Shift right by 8 bits to get the sign bit
         
         # Apply sign to current values
+        # We need to cast to int16 to handle negative values correctly
+        # The max value of 8 bits is 255 so the casting is safe as there are non-negative values
+        magnitude = magnitude.astype(np.int16)
         current = np.where(sign_bit == 1, -magnitude, magnitude)
         
         # Note: If needed, other flag bits could be extracted as follows:

From b5075629f58711777e10f121a4d6fb6d30abb300 Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Fri, 21 Mar 2025 12:10:14 -0600
Subject: [PATCH 08/20] refactor and comments

---
 neo/rawio/intanrawio.py | 76 ++++++++++++++++++++---------------------
 1 file changed, 37 insertions(+), 39 deletions(-)

diff --git a/neo/rawio/intanrawio.py b/neo/rawio/intanrawio.py
index 4defea300..0289269af 100644
--- a/neo/rawio/intanrawio.py
+++ b/neo/rawio/intanrawio.py
@@ -462,28 +462,22 @@ def _get_analogsignal_chunk_header_attached(self, i_start, i_stop, stream_index,
             sl0 = i_start % block_size
             sl1 = sl0 + (i_stop - i_start)
 
-        # For all streams raw_data is a structured memmap with a field for each channel_id
-        if stream_is_stim:
-            # For stim data, we need to extract the raw data first, then demultiplex it
-            stim_data = np.zeros((i_stop - i_start, len(channel_ids)), dtype=dtype)
-            for chunk_index, channel_id in enumerate(channel_ids):
-                data_chan = self._raw_data[channel_id]
-                if multiple_samples_per_block:
-                    stim_data[:, chunk_index] = data_chan[block_start:block_stop].flatten()[sl0:sl1]
-                else:
-                    stim_data[:, chunk_index] = data_chan[i_start:i_stop]
-            # Now demultiplex the stim data
-            sigs_chunk = self._decode_current_from_stim_data(stim_data, 0, stim_data.shape[0])
-        elif not stream_is_digital:
-            sigs_chunk = np.zeros((i_stop - i_start, len(channel_ids)), dtype=dtype)
-
+        if not stream_is_digital:
+            # For all streams raw_data is a structured memmap with a field for each channel_id        
+            sigs_chunk = np.zeros((i_stop - i_start, len(channel_ids)), dtype=dtype)            
             for chunk_index, channel_id in enumerate(channel_ids):
                 data_chan = self._raw_data[channel_id]
+                
                 if multiple_samples_per_block:
                     sigs_chunk[:, chunk_index] = data_chan[block_start:block_stop].flatten()[sl0:sl1]
                 else:
                     sigs_chunk[:, chunk_index] = data_chan[i_start:i_stop]
-        else:  # For digital data the channels come interleaved in a single field so we need to demultiplex
+
+            if stream_is_stim:
+                sigs_chunk = self._decode_current_from_stim_data(sigs_chunk, 0, sigs_chunk.shape[0])
+
+        else: 
+             # For digital data the channels come interleaved in a single field so we need to demultiplex
             digital_raw_data = self._raw_data[field_name].flatten()
             sigs_chunk = self._demultiplex_digital_data(digital_raw_data, channel_ids, i_start, i_stop)
         return sigs_chunk
@@ -493,7 +487,7 @@ def _get_analogsignal_chunk_one_file_per_channel(self, i_start, i_stop, stream_i
         stream_name = self.header["signal_streams"][stream_index]["name"][:]
         signal_data_memmap_list = self._raw_data[stream_name]
         stream_is_stim = stream_name == "Stim channel"
-        
+
         channel_indexes_are_slice = isinstance(channel_indexes, slice)
         if channel_indexes_are_slice:
             num_channels = len(signal_data_memmap_list)
@@ -510,8 +504,8 @@ def _get_analogsignal_chunk_one_file_per_channel(self, i_start, i_stop, stream_i
         for chunk_index, channel_index in enumerate(channel_indexes):
             channel_memmap = signal_data_memmap_list[channel_index]
             sigs_chunk[:, chunk_index] = channel_memmap[i_start:i_stop]
-            
-        # If this is stim data, we need to demultiplex it
+
+        # If this is stim data, we need to decode the current values
         if stream_is_stim:
             sigs_chunk = self._decode_current_from_stim_data(sigs_chunk, 0, sigs_chunk.shape[0])
 
@@ -524,7 +518,7 @@ def _get_analogsignal_chunk_one_file_per_signal(self, i_start, i_stop, stream_in
 
         stream_is_digital = stream_name in digital_stream_names
         stream_is_stim = stream_name == "Stim channel"
-        
+
         if stream_is_digital:
             stream_id = self.header["signal_streams"][stream_index]["id"]
             mask = self.header["signal_channels"]["stream_id"] == stream_id
@@ -543,11 +537,11 @@ def _get_analogsignal_chunk_one_file_per_signal(self, i_start, i_stop, stream_in
     def _demultiplex_digital_data(self, raw_digital_data, channel_ids, i_start, i_stop):
         """
         Demultiplex digital data by extracting individual channel values from packed 16-bit format.
-        
+
         According to the Intan format, digital input/output data is stored with all 16 channels
         encoded bit-by-bit in each 16-bit word. This method extracts the specified digital channels
         from the packed format into separate uint16 arrays of 0 and 1.
-        
+
         Parameters
         ----------
         raw_digital_data : ndarray
@@ -559,19 +553,19 @@ def _demultiplex_digital_data(self, raw_digital_data, channel_ids, i_start, i_st
             Starting sample index for demultiplexing.
         i_stop : int
             Ending sample index for demultiplexing (exclusive).
-            
+
         Returns
         -------
         ndarray
-            Demultiplexed digital data with shape (i_stop-i_start, len(channel_ids)), 
+            Demultiplexed digital data with shape (i_stop-i_start, len(channel_ids)),
             containing 0 or 1 values for each requested channel.
-        
+
         Notes
         -----
         In the Intan format, digital channels are packed into 16-bit words where each bit position
         corresponds to a specific channel number. For example, with digital inputs 0, 4, and 5
         set high and the rest low, the 16-bit word would be 2^0 + 2^4 + 2^5 = 1 + 16 + 32 = 49.
-        
+
         The native_order property for each channel corresponds to its bit position in the packed word.
         """
         dtype = np.uint16  # We fix this to match the memmap dtype
@@ -584,15 +578,15 @@ def _demultiplex_digital_data(self, raw_digital_data, channel_ids, i_start, i_st
             output[:, channel_index] = demultiplex_data[i_start:i_stop].flatten()
 
         return output
-        
+
     def _decode_current_from_stim_data(self, raw_stim_data, i_start, i_stop):
         """
         Demultiplex stimulation data by extracting current values from packed 16-bit format.
-        
+
         According to the Intan RHS data format, stimulation current is stored in the lower 9 bits
         of each 16-bit word: 8 bits for magnitude and 1 bit for sign. The upper bits contain
         flags for compliance limit, charge recovery, and amplifier settle.
-        
+
         Parameters
         ----------
         raw_stim_data : ndarray
@@ -601,14 +595,14 @@ def _decode_current_from_stim_data(self, raw_stim_data, i_start, i_stop):
             Starting sample index for demultiplexing.
         i_stop : int
             Ending sample index for demultiplexing (exclusive).
-            
+
         Returns
         -------
         ndarray
             Demultiplexed stimulation current values in amperes, preserving the original
             array dimensions. The output values need to be multiplied by the stim_step_size
             parameter (from header) to obtain the actual current in amperes.
-        
+
         Notes
         -----
         Bit structure of each 16-bit stimulation word:
@@ -616,32 +610,32 @@ def _decode_current_from_stim_data(self, raw_stim_data, i_start, i_stop):
         - Bit 8: Sign bit (1 = negative current)
         - Bits 9-13: Unused (always zero)
         - Bit 14: Amplifier settle flag (1 = activated)
-        - Bit 15: Charge recovery flag (1 = activated) 
+        - Bit 15: Charge recovery flag (1 = activated)
         - Bit 16 (MSB): Compliance limit flag (1 = limit reached)
-        
+
         The actual current value in amperes is obtained by multiplying the
-        output by the 'stim_step_size' parameter from the file header. These scaled values can be 
+        output by the 'stim_step_size' parameter from the file header. These scaled values can be
         obtained with the `rescale_signal_raw_to_float` function.
         """
         # Get the relevant portion of the data
         data = raw_stim_data[i_start:i_stop]
-        
+
         # Extract current value (bits 0-8)
         magnitude = np.bitwise_and(data, 0xFF)
         sign_bit = np.bitwise_and(np.right_shift(data, 8), 0x01)  # Shift right by 8 bits to get the sign bit
-        
+
         # Apply sign to current values
         # We need to cast to int16 to handle negative values correctly
         # The max value of 8 bits is 255 so the casting is safe as there are non-negative values
         magnitude = magnitude.astype(np.int16)
         current = np.where(sign_bit == 1, -magnitude, magnitude)
-        
+
         # Note: If needed, other flag bits could be extracted as follows:
         # compliance_flag = np.bitwise_and(np.right_shift(data, 15), 0x01).astype(bool)  # Bit 16 (MSB)
         # charge_recovery_flag = np.bitwise_and(np.right_shift(data, 14), 0x01).astype(bool)  # Bit 15
         # amp_settle_flag = np.bitwise_and(np.right_shift(data, 13), 0x01).astype(bool)  # Bit 14
         # These could be returned as a structured array or dictionary if needed
-        
+
         return current
 
     def get_intan_timestamps(self, i_start=None, i_stop=None):
@@ -974,8 +968,12 @@ def read_rhs(filename, file_format: str):
             chan_info_stim["gain"] = global_info["stim_step_size"]
             chan_info_stim["offset"] = 0.0
             chan_info_stim["signal_type"] = 11  # put it in another group
-            chan_info_stim["dtype"] = "uint16"
+            chan_info_stim["dtype"] = "int16"
             ordered_channel_info.append(chan_info_stim)
+            
+            # Note that the data on disk is uint16 but the data is 
+            # then decoded as int16 so the chan_info is int16
+            memmap_dtype = "uint16"
             if file_format == "header-attached":
                 memmap_data_dtype += [(name + "_STIM", "uint16", BLOCK_SIZE)]
             else:

From 17138ddaeb3ac533f380e1262a888a50ed8f94ec Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Mon, 24 Mar 2025 12:56:25 -0600
Subject: [PATCH 09/20] add test

---
 neo/test/iotest/test_intanio.py       |  1 +
 neo/test/rawiotest/test_intanrawio.py | 50 +++++++++++++++++++++++++++
 2 files changed, 51 insertions(+)

diff --git a/neo/test/iotest/test_intanio.py b/neo/test/iotest/test_intanio.py
index 1e527f052..ccd3be90c 100644
--- a/neo/test/iotest/test_intanio.py
+++ b/neo/test/iotest/test_intanio.py
@@ -23,6 +23,7 @@ class TestIntanIO(
         "intan/intan_fpc_rhs_test_240329_091637/info.rhs",  # Format one-file-per-channel
         "intan/intan_fps_rhs_test_240329_091536/info.rhs",  # Format one-file-per-signal
         "intan/rhd_fpc_multistim_240514_082044/info.rhd",  # Multiple digital channels one-file-per-channel rhd
+        "intan/rhs_stim_data_single_file_format/intanTestFile.rhs", # header-attached rhs data with stimulus current 
     ]
 
 
diff --git a/neo/test/rawiotest/test_intanrawio.py b/neo/test/rawiotest/test_intanrawio.py
index 3d86e5089..ab90ec586 100644
--- a/neo/test/rawiotest/test_intanrawio.py
+++ b/neo/test/rawiotest/test_intanrawio.py
@@ -21,6 +21,7 @@ class TestIntanRawIO(
         "intan/intan_fpc_rhs_test_240329_091637/info.rhs",  # Format one-file-per-channel
         "intan/intan_fps_rhs_test_240329_091536/info.rhs",  # Format one-file-per-signal
         "intan/rhd_fpc_multistim_240514_082044/info.rhd",  # Multiple digital channels one-file-per-channel rhd
+        "intan/rhs_stim_data_single_file_format/intanTestFile.rhs", # header-attached rhs data with stimulus current 
     ]
 
     def test_annotations(self):
@@ -87,5 +88,54 @@ def test_correct_reading_one_file_per_channel(self):
             np.testing.assert_allclose(data_raw, data_from_neo)
 
 
+    def test_correct_decoding_of_stimulus_current(self):
+        
+        file_path = Path(self.get_local_path("intan/rhs_stim_data_single_file_format/intanTestFile.rhs"))
+        intan_reader = IntanRawIO(filename=file_path)
+        intan_reader.parse_header()
+        
+        signal_streams = intan_reader.header['signal_streams']
+        stream_ids = signal_streams['id'].tolist()
+        stream_index = stream_ids.index('11')
+        sampling_rate = intan_reader.get_signal_sampling_rate(stream_index=stream_index)
+        sig_chunk = intan_reader.get_analogsignal_chunk(stream_index=stream_index, channel_ids=["D-016_STIM"])
+        final_stim = intan_reader.rescale_signal_raw_to_float(sig_chunk, stream_index=stream_index, channel_ids=["D-016_STIM"])
+
+        # This contains only the first pulse and I got this by visual inspection
+        data_to_test = final_stim[200:250]
+        
+        positive_pulse_size = np.max(data_to_test).item()
+        negative_pulse_size = np.min(data_to_test).item()
+
+        expected_value = 60 * 1e-6# 60 microamperes
+
+        # Assert is close float
+        assert np.isclose(positive_pulse_size, expected_value)
+        assert np.isclose(negative_pulse_size, -expected_value)
+
+        # Check that negative pulse is leading
+        argmin = np.argmin(data_to_test)
+        argmax = np.argmax(data_to_test)
+
+        assert argmin < argmax
+
+
+        # Check that the negative pulse is 200 us long
+        negative_pulse_frames = np.where(data_to_test > 0)[0]
+        number_of_negative_frames = negative_pulse_frames.size
+        duration_of_negative_pulse = number_of_negative_frames / sampling_rate
+
+        expected_duration = 200 * 1e-6  # 400 microseconds / 2 
+
+        assert np.isclose(duration_of_negative_pulse, expected_duration, rtol=1e-05, atol=1e-08)
+
+        positive_pulse_frames = np.where(data_to_test > 0)[0]
+        number_of_positive_frames = positive_pulse_frames.size
+        duration_of_positive_pulse = number_of_positive_frames / sampling_rate
+        expected_duration = 200 * 1e-6  # 400 microseconds / 2
+
+        assert np.isclose(duration_of_positive_pulse, expected_duration, rtol=1e-05, atol=1e-08)
+
+
 if __name__ == "__main__":
     unittest.main()

From af6f62881af0a76742d7037a115eaef24beb3e20 Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Mon, 24 Mar 2025 13:02:22 -0600
Subject: [PATCH 10/20] add comments to test

---
 neo/test/rawiotest/test_intanrawio.py | 13 ++++++-------
 1 file changed, 6 insertions(+), 7 deletions(-)

diff --git a/neo/test/rawiotest/test_intanrawio.py b/neo/test/rawiotest/test_intanrawio.py
index ab90ec586..9911d7431 100644
--- a/neo/test/rawiotest/test_intanrawio.py
+++ b/neo/test/rawiotest/test_intanrawio.py
@@ -89,6 +89,7 @@ def test_correct_reading_one_file_per_channel(self):
 
 
     def test_correct_decoding_of_stimulus_current(self):
+        # See https://github.com/NeuralEnsemble/python-neo/pull/1660 for discussion
         
         file_path = Path(self.get_local_path("intan/rhs_stim_data_single_file_format/intanTestFile.rhs"))
         intan_reader = IntanRawIO(filename=file_path)
@@ -116,25 +117,23 @@ def test_correct_decoding_of_stimulus_current(self):
         # Check that negative pulse is leading
         argmin = np.argmin(data_to_test)
         argmax = np.argmax(data_to_test)
-
         assert argmin < argmax
-
-
+        
         # Check that the negative pulse is 200 us long
         negative_pulse_frames = np.where(data_to_test > 0)[0]
         number_of_negative_frames = negative_pulse_frames.size
         duration_of_negative_pulse = number_of_negative_frames / sampling_rate
 
         expected_duration = 200 * 1e-6  # 400 microseconds / 2 
-
-        assert np.isclose(duration_of_negative_pulse, expected_duration, rtol=1e-05, atol=1e-08)
-
+        assert np.isclose(duration_of_negative_pulse, expected_duration)
+        
+        # Check that the positive pulse is 200 us long
         positive_pulse_frames = np.where(data_to_test > 0)[0]
         number_of_positive_frames = positive_pulse_frames.size
         duration_of_positive_pulse = number_of_positive_frames / sampling_rate
         expected_duration = 200 * 1e-6  # 400 microseconds / 2
 
-        assert np.isclose(duration_of_positive_pulse, expected_duration, rtol=1e-05, atol=1e-08)
+        assert np.isclose(duration_of_positive_pulse, expected_duration)
 
 
 if __name__ == "__main__":

From 742ebca55298d178f4a240be8c376b5acc8f98f7 Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Mon, 24 Mar 2025 13:03:15 -0600
Subject: [PATCH 11/20] add gin comment

---
 neo/test/rawiotest/test_intanrawio.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/neo/test/rawiotest/test_intanrawio.py b/neo/test/rawiotest/test_intanrawio.py
index 9911d7431..c4a355363 100644
--- a/neo/test/rawiotest/test_intanrawio.py
+++ b/neo/test/rawiotest/test_intanrawio.py
@@ -90,6 +90,8 @@ def test_correct_reading_one_file_per_channel(self):
 
     def test_correct_decoding_of_stimulus_current(self):
         # See https://github.com/NeuralEnsemble/python-neo/pull/1660 for discussion
+        # See https://gin.g-node.org/NeuralEnsemble/ephy_testing_data/src/master/intan/README.md#rhs_stim_data_single_file_format 
+        # For a description of the data 
         
         file_path = Path(self.get_local_path("intan/rhs_stim_data_single_file_format/intanTestFile.rhs"))
         intan_reader = IntanRawIO(filename=file_path)

From 1eac371408d06042d1a657bc7ea59faa47f1f35e Mon Sep 17 00:00:00 2001
From: zm711 <92116279+zm711@users.noreply.github.com>
Date: Thu, 27 Mar 2025 14:18:43 -0400
Subject: [PATCH 12/20] update header parsing for stim + update docstring

---
 neo/rawio/intanrawio.py | 70 +++++++++++++++++++++--------------------
 1 file changed, 36 insertions(+), 34 deletions(-)

diff --git a/neo/rawio/intanrawio.py b/neo/rawio/intanrawio.py
index 4cda50071..31d2cdc4e 100644
--- a/neo/rawio/intanrawio.py
+++ b/neo/rawio/intanrawio.py
@@ -63,14 +63,12 @@ class IntanRawIO(BaseRawIO):
       'one-file-per-channel' which have a header file called 'info.rhd' or 'info.rhs' and a series
       of binary files with the '.dat' suffix
 
-    * The reader can handle three file formats 'header-attached', 'one-file-per-signal' and
-      'one-file-per-channel'.
-
-    * Intan files contain amplifier channels labeled 'A', 'B' 'C' or 'D'
-      depending on the port in which they were recorded along with the following
+    * Intan files contain amplifier channels labeled 'A', 'B' 'C' or 'D' for the 512 recorder
+      or 'A', 'B', 'C', 'D', 'E', 'F', 'G', 'H' for the 1024 recorder system
+      depending on the port in which they were recorded along (stored in stream_id '0') with the following
       additional streams.
 
-    0: 'RHD2000' amplifier channel
+    0: 'RHD2000 amplifier channel'
     1: 'RHD2000 auxiliary input channel',
     2: 'RHD2000 supply voltage channel',
     3: 'USB board ADC input channel',
@@ -87,9 +85,11 @@ class IntanRawIO(BaseRawIO):
     10: 'DC Amplifier channel',
     11: 'Stim channel',
 
-    * For the "header-attached" and "one-file-per-signal" formats, the structure of the digital input and output channels is
-      one long vector, which must be post-processed to extract individual digital channel information.
-      See the intantech website for more information on performing this post-processing.
+    * We currently implement digital data demultiplexing so that if digital streams are requested they are
+      returned as arrays of 1s and 0s.
+
+    * We also do stim data decoding which returns the stim data as an int16 of appropriate magnitude. Please
+      use `rescale_signal_raw_to_float` to obtain stim data in amperes.
 
 
     Examples
@@ -954,31 +954,33 @@ def read_rhs(filename, file_format: str):
                 memmap_data_dtype["DC Amplifier channel"] = "uint16"
 
     # I can't seem to get stim files to generate for one-file-per-channel
-    # so let's skip for now and can be given on request
-    if file_format != "one-file-per-channel":
-        for chan_info in stream_name_to_channel_info_list["RHS2000 amplifier channel"]:
-            chan_info_stim = dict(chan_info)
-            name = chan_info["native_channel_name"]
-            chan_info_stim["native_channel_name"] = name + "_STIM"
-            chan_info_stim["sampling_rate"] = sr
-            # stim channel are complicated because they are coded
-            # with bits, they do not fit the gain/offset rawio strategy
-            chan_info_stim["units"] = "A"  # Amps
-            chan_info_stim["gain"] = global_info["stim_step_size"]
-            chan_info_stim["offset"] = 0.0
-            chan_info_stim["signal_type"] = 11  # put it in another group
-            chan_info_stim["dtype"] = "int16"
-            ordered_channel_info.append(chan_info_stim)
-            
-            # Note that the data on disk is uint16 but the data is 
-            # then decoded as int16 so the chan_info is int16
-            memmap_dtype = "uint16"
-            if file_format == "header-attached":
-                memmap_data_dtype += [(name + "_STIM", "uint16", BLOCK_SIZE)]
-            else:
-                memmap_data_dtype["Stim channel"] = "uint16"
-    else:
-        warnings.warn("Stim not implemented for `one-file-per-channel` due to lack of test files")
+    # so ideally at some point we need test data to confirm this is true
+    # based on what Heberto and I read in the docs
+    for chan_info in stream_name_to_channel_info_list["RHS2000 amplifier channel"]:
+        chan_info_stim = dict(chan_info)
+        name = chan_info["native_channel_name"]
+        chan_info_stim["native_channel_name"] = name + "_STIM"
+        chan_info_stim["sampling_rate"] = sr
+        # stim channel are complicated because they are coded
+        # with bits, they do not fit the gain/offset rawio strategy
+        chan_info_stim["units"] = "A"  # Amps
+        chan_info_stim["gain"] = global_info["stim_step_size"]
+        chan_info_stim["offset"] = 0.0
+        chan_info_stim["signal_type"] = 11  # put it in another group
+        chan_info_stim["dtype"] = "int16" # this change is due to bit decoding see note below
+        ordered_channel_info.append(chan_info_stim)
+        # Note that the data on disk is uint16 but the data is 
+        # then decoded as int16 so the chan_info is int16
+        if file_format == "header-attached":
+            memmap_data_dtype += [(name + "_STIM", "uint16", BLOCK_SIZE)]
+        else:
+            memmap_data_dtype["Stim channel"] = "uint16"
+        if file_format == "one-file-per-channel":
+            warning_msg = ("Stim decoding for `one-file-per-channel` is based on a reading of the documention "
+                           "and we would appreciate test data to ensure we've implemented this format "
+                           "appropriately. If you use the stim data please verify it is as you expected "
+                           "and if it is not then open an issue on the python-neo repo")
+            warnings.warn(warning_msg)
 
     # No supply or aux for rhs files (ie no stream_id 1 and 2)
     # We have an error above that requests test files to help if the spec is changed

From e45e1fc4f69acbb15ac00642cbaf9bf62ca62f1a Mon Sep 17 00:00:00 2001
From: zm711 <92116279+zm711@users.noreply.github.com>
Date: Thu, 27 Mar 2025 14:32:11 -0400
Subject: [PATCH 13/20] add stim to channel check

---
 neo/rawio/intanrawio.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/neo/rawio/intanrawio.py b/neo/rawio/intanrawio.py
index 31d2cdc4e..047cc2ebf 100644
--- a/neo/rawio/intanrawio.py
+++ b/neo/rawio/intanrawio.py
@@ -1446,6 +1446,7 @@ def create_one_file_per_channel_dict_rhd(dirname):
     "USB board ADC output channel": "board-ANALOG-OUT",
     "USB board digital input channel": "board-DIGITAL-IN",
     "USB board digital output channel": "board-DIGITAL-OUT",
+    'Stim channel': "stim"
 }
 
 

From d5737271ed91110457ec1b4b4ccd7dfc215829fc Mon Sep 17 00:00:00 2001
From: zm711 <92116279+zm711@users.noreply.github.com>
Date: Thu, 27 Mar 2025 16:51:11 -0400
Subject: [PATCH 14/20] wip

---
 neo/rawio/intanrawio.py | 92 ++++++++++++++++++++++-------------------
 1 file changed, 49 insertions(+), 43 deletions(-)

diff --git a/neo/rawio/intanrawio.py b/neo/rawio/intanrawio.py
index 047cc2ebf..752e834d0 100644
--- a/neo/rawio/intanrawio.py
+++ b/neo/rawio/intanrawio.py
@@ -900,8 +900,17 @@ def read_rhs(filename, file_format: str):
         channel_number_dict = {name: len(stream_name_to_channel_info_list[name]) for name in names_to_count}
 
         # Both DC Amplifier and Stim streams have the same number of channels as the amplifier stream
+        # if using the `header-attached` or `one-file-per-signal` formats
+        # the amplifier data is stored in the same place in the header so no matter what the DC amp
+        # uses the same info as the RHS amp.
         channel_number_dict["DC Amplifier channel"] = channel_number_dict["RHS2000 amplifier channel"]
-        channel_number_dict["Stim channel"] = channel_number_dict["RHS2000 amplifier channel"]
+        if file_format != "one-file-per-channel":
+            channel_number_dict["Stim channel"] = channel_number_dict["RHS2000 amplifier channel"]
+        else:
+            raw_file_paths_dict = create_one_file_per_channel_dict_rhs(dirname=filename.parent)
+            channel_number_dict["Stim channel"] = len(raw_file_paths_dict["Stim channel"])
+            # but the user can shut off the normal amplifier and only save dc amplifier
+            channel_number_dict["RHS2000 amplifier channel"] = len(raw_file_paths_dict["RHS2000 amplifier channel"])
 
         header_size = f.tell()
 
@@ -915,24 +924,25 @@ def read_rhs(filename, file_format: str):
         memmap_data_dtype["timestamp"] = "int32"
         channel_number_dict["timestamp"] = 1
 
-    for chan_info in stream_name_to_channel_info_list["RHS2000 amplifier channel"]:
-        chan_info["sampling_rate"] = sr
-        chan_info["units"] = "uV"
-        chan_info["gain"] = 0.195
-        if file_format == "header-attached":
-            chan_info["offset"] = -32768 * 0.195
-        else:
-            chan_info["offset"] = 0.0
-        if file_format == "header-attached":
-            chan_info["dtype"] = "uint16"
-        else:
-            chan_info["dtype"] = "int16"
-        ordered_channel_info.append(chan_info)
-        if file_format == "header-attached":
-            name = chan_info["native_channel_name"]
-            memmap_data_dtype += [(name, "uint16", BLOCK_SIZE)]
-        else:
-            memmap_data_dtype["RHS2000 amplifier channel"] = "int16"
+    if file_format != "one-file-per-channel" or channel_number_dict["RHS2000 amplifier channel"] > 0:
+        for chan_info in stream_name_to_channel_info_list["RHS2000 amplifier channel"]:
+            chan_info["sampling_rate"] = sr
+            chan_info["units"] = "uV"
+            chan_info["gain"] = 0.195
+            if file_format == "header-attached":
+                chan_info["offset"] = -32768 * 0.195
+            else:
+                chan_info["offset"] = 0.0
+            if file_format == "header-attached":
+                chan_info["dtype"] = "uint16"
+            else:
+                chan_info["dtype"] = "int16"
+            ordered_channel_info.append(chan_info)
+            if file_format == "header-attached":
+                name = chan_info["native_channel_name"]
+                memmap_data_dtype += [(name, "uint16", BLOCK_SIZE)]
+            else:
+                memmap_data_dtype["RHS2000 amplifier channel"] = "int16"
 
     if bool(global_info["dc_amplifier_data_saved"]):
         # if we have dc amp we need to grab the correct number of channels
@@ -957,30 +967,26 @@ def read_rhs(filename, file_format: str):
     # so ideally at some point we need test data to confirm this is true
     # based on what Heberto and I read in the docs
     for chan_info in stream_name_to_channel_info_list["RHS2000 amplifier channel"]:
-        chan_info_stim = dict(chan_info)
-        name = chan_info["native_channel_name"]
-        chan_info_stim["native_channel_name"] = name + "_STIM"
-        chan_info_stim["sampling_rate"] = sr
-        # stim channel are complicated because they are coded
-        # with bits, they do not fit the gain/offset rawio strategy
-        chan_info_stim["units"] = "A"  # Amps
-        chan_info_stim["gain"] = global_info["stim_step_size"]
-        chan_info_stim["offset"] = 0.0
-        chan_info_stim["signal_type"] = 11  # put it in another group
-        chan_info_stim["dtype"] = "int16" # this change is due to bit decoding see note below
-        ordered_channel_info.append(chan_info_stim)
-        # Note that the data on disk is uint16 but the data is 
-        # then decoded as int16 so the chan_info is int16
-        if file_format == "header-attached":
-            memmap_data_dtype += [(name + "_STIM", "uint16", BLOCK_SIZE)]
-        else:
-            memmap_data_dtype["Stim channel"] = "uint16"
-        if file_format == "one-file-per-channel":
-            warning_msg = ("Stim decoding for `one-file-per-channel` is based on a reading of the documention "
-                           "and we would appreciate test data to ensure we've implemented this format "
-                           "appropriately. If you use the stim data please verify it is as you expected "
-                           "and if it is not then open an issue on the python-neo repo")
-            warnings.warn(warning_msg)
+        # we see which stim were activated
+        if any([chan_info["native_channel_name"] in stim_file.stem for stim_file in raw_file_paths_dict['Stim channel']]):
+            chan_info_stim = dict(chan_info)
+            name = chan_info["native_channel_name"]
+            chan_info_stim["native_channel_name"] = name + "_STIM"
+            chan_info_stim["sampling_rate"] = sr
+            # stim channel are complicated because they are coded
+            # with bits, they do not fit the gain/offset rawio strategy
+            chan_info_stim["units"] = "A"  # Amps
+            chan_info_stim["gain"] = global_info["stim_step_size"]
+            chan_info_stim["offset"] = 0.0
+            chan_info_stim["signal_type"] = 11  # put it in another group
+            chan_info_stim["dtype"] = "int16" # this change is due to bit decoding see note below
+            ordered_channel_info.append(chan_info_stim)
+            # Note that the data on disk is uint16 but the data is 
+            # then decoded as int16 so the chan_info is int16
+            if file_format == "header-attached":
+                memmap_data_dtype += [(name + "_STIM", "uint16", BLOCK_SIZE)]
+            else:
+                memmap_data_dtype["Stim channel"] = "uint16"
 
     # No supply or aux for rhs files (ie no stream_id 1 and 2)
     # We have an error above that requests test files to help if the spec is changed

From a83f73cc319387c5ebf4af7f2bc431da588cb121 Mon Sep 17 00:00:00 2001
From: zm711 <92116279+zm711@users.noreply.github.com>
Date: Thu, 27 Mar 2025 16:57:16 -0400
Subject: [PATCH 15/20] oops

---
 neo/rawio/intanrawio.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/neo/rawio/intanrawio.py b/neo/rawio/intanrawio.py
index 752e834d0..c3573d1fc 100644
--- a/neo/rawio/intanrawio.py
+++ b/neo/rawio/intanrawio.py
@@ -906,6 +906,7 @@ def read_rhs(filename, file_format: str):
         channel_number_dict["DC Amplifier channel"] = channel_number_dict["RHS2000 amplifier channel"]
         if file_format != "one-file-per-channel":
             channel_number_dict["Stim channel"] = channel_number_dict["RHS2000 amplifier channel"]
+            raw_file_paths_dict = create_one_file_per_signal_dict_rhs(dirname=filename.parent)
         else:
             raw_file_paths_dict = create_one_file_per_channel_dict_rhs(dirname=filename.parent)
             channel_number_dict["Stim channel"] = len(raw_file_paths_dict["Stim channel"])

From 721e11935371b3f9e785ccdf244483af12eb8ff7 Mon Sep 17 00:00:00 2001
From: zm711 <92116279+zm711@users.noreply.github.com>
Date: Thu, 27 Mar 2025 17:13:38 -0400
Subject: [PATCH 16/20] fix for header attached

---
 neo/rawio/intanrawio.py | 19 ++++++++++---------
 1 file changed, 10 insertions(+), 9 deletions(-)

diff --git a/neo/rawio/intanrawio.py b/neo/rawio/intanrawio.py
index c3573d1fc..22f589ea5 100644
--- a/neo/rawio/intanrawio.py
+++ b/neo/rawio/intanrawio.py
@@ -462,11 +462,11 @@ def _get_analogsignal_chunk_header_attached(self, i_start, i_stop, stream_index,
             sl1 = sl0 + (i_stop - i_start)
 
         if not stream_is_digital:
-            # For all streams raw_data is a structured memmap with a field for each channel_id        
-            sigs_chunk = np.zeros((i_stop - i_start, len(channel_ids)), dtype=dtype)            
+            # For all streams raw_data is a structured memmap with a field for each channel_id
+            sigs_chunk = np.zeros((i_stop - i_start, len(channel_ids)), dtype=dtype)
             for chunk_index, channel_id in enumerate(channel_ids):
                 data_chan = self._raw_data[channel_id]
-                
+
                 if multiple_samples_per_block:
                     sigs_chunk[:, chunk_index] = data_chan[block_start:block_stop].flatten()[sl0:sl1]
                 else:
@@ -475,8 +475,8 @@ def _get_analogsignal_chunk_header_attached(self, i_start, i_stop, stream_index,
             if stream_is_stim:
                 sigs_chunk = self._decode_current_from_stim_data(sigs_chunk, 0, sigs_chunk.shape[0])
 
-        else: 
-             # For digital data the channels come interleaved in a single field so we need to demultiplex
+        else:
+            # For digital data the channels come interleaved in a single field so we need to demultiplex
             digital_raw_data = self._raw_data[field_name].flatten()
             sigs_chunk = self._demultiplex_digital_data(digital_raw_data, channel_ids, i_start, i_stop)
         return sigs_chunk
@@ -969,7 +969,9 @@ def read_rhs(filename, file_format: str):
     # based on what Heberto and I read in the docs
     for chan_info in stream_name_to_channel_info_list["RHS2000 amplifier channel"]:
         # we see which stim were activated
-        if any([chan_info["native_channel_name"] in stim_file.stem for stim_file in raw_file_paths_dict['Stim channel']]):
+        if file_format == "header-attached" or any(
+            [chan_info["native_channel_name"] in stim_file.stem for stim_file in raw_file_paths_dict["Stim channel"]]
+        ):
             chan_info_stim = dict(chan_info)
             name = chan_info["native_channel_name"]
             chan_info_stim["native_channel_name"] = name + "_STIM"
@@ -980,9 +982,9 @@ def read_rhs(filename, file_format: str):
             chan_info_stim["gain"] = global_info["stim_step_size"]
             chan_info_stim["offset"] = 0.0
             chan_info_stim["signal_type"] = 11  # put it in another group
-            chan_info_stim["dtype"] = "int16" # this change is due to bit decoding see note below
+            chan_info_stim["dtype"] = "int16"  # this change is due to bit decoding see note below
             ordered_channel_info.append(chan_info_stim)
-            # Note that the data on disk is uint16 but the data is 
+            # Note that the data on disk is uint16 but the data is
             # then decoded as int16 so the chan_info is int16
             if file_format == "header-attached":
                 memmap_data_dtype += [(name + "_STIM", "uint16", BLOCK_SIZE)]
@@ -1453,7 +1455,6 @@ def create_one_file_per_channel_dict_rhd(dirname):
     "USB board ADC output channel": "board-ANALOG-OUT",
     "USB board digital input channel": "board-DIGITAL-IN",
     "USB board digital output channel": "board-DIGITAL-OUT",
-    'Stim channel': "stim"
 }
 
 

From e05d7722d7efdc5ae1dd61cc4e371e6b26e3b8dd Mon Sep 17 00:00:00 2001
From: zm711 <92116279+zm711@users.noreply.github.com>
Date: Thu, 27 Mar 2025 17:20:39 -0400
Subject: [PATCH 17/20] fix for one-file-per-signal

---
 neo/rawio/intanrawio.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/neo/rawio/intanrawio.py b/neo/rawio/intanrawio.py
index 22f589ea5..6bd466d43 100644
--- a/neo/rawio/intanrawio.py
+++ b/neo/rawio/intanrawio.py
@@ -969,7 +969,7 @@ def read_rhs(filename, file_format: str):
     # based on what Heberto and I read in the docs
     for chan_info in stream_name_to_channel_info_list["RHS2000 amplifier channel"]:
         # we see which stim were activated
-        if file_format == "header-attached" or any(
+        if file_format != "one-file-per-channel" or any(
             [chan_info["native_channel_name"] in stim_file.stem for stim_file in raw_file_paths_dict["Stim channel"]]
         ):
             chan_info_stim = dict(chan_info)

From 42d56fdf7ef660e5fa25bcb314dd871e557afa9e Mon Sep 17 00:00:00 2001
From: zm711 <92116279+zm711@users.noreply.github.com>
Date: Sat, 29 Mar 2025 08:53:23 -0400
Subject: [PATCH 18/20] add one new test file

---
 neo/test/iotest/test_intanio.py       | 2 ++
 neo/test/rawiotest/test_intanrawio.py | 4 +++-
 2 files changed, 5 insertions(+), 1 deletion(-)

diff --git a/neo/test/iotest/test_intanio.py b/neo/test/iotest/test_intanio.py
index ccd3be90c..f76659338 100644
--- a/neo/test/iotest/test_intanio.py
+++ b/neo/test/iotest/test_intanio.py
@@ -24,6 +24,8 @@ class TestIntanIO(
         "intan/intan_fps_rhs_test_240329_091536/info.rhs",  # Format one-file-per-signal
         "intan/rhd_fpc_multistim_240514_082044/info.rhd",  # Multiple digital channels one-file-per-channel rhd
         "intan/rhs_stim_data_single_file_format/intanTestFile.rhs", # header-attached rhs data with stimulus current 
+        "intan/test_fcs_dc_250327_154333/info.rhs", # this is an example of only having dc amp rather than amp files
+        #"intan/test_fpc_stim_250327_151617/info.rhs", # wrong files Heberto will fix
     ]
 
 
diff --git a/neo/test/rawiotest/test_intanrawio.py b/neo/test/rawiotest/test_intanrawio.py
index c4a355363..1ccf7b911 100644
--- a/neo/test/rawiotest/test_intanrawio.py
+++ b/neo/test/rawiotest/test_intanrawio.py
@@ -21,7 +21,9 @@ class TestIntanRawIO(
         "intan/intan_fpc_rhs_test_240329_091637/info.rhs",  # Format one-file-per-channel
         "intan/intan_fps_rhs_test_240329_091536/info.rhs",  # Format one-file-per-signal
         "intan/rhd_fpc_multistim_240514_082044/info.rhd",  # Multiple digital channels one-file-per-channel rhd
-        "intan/rhs_stim_data_single_file_format/intanTestFile.rhs", # header-attached rhs data with stimulus current 
+        "intan/rhs_stim_data_single_file_format/intanTestFile.rhs", # header-attached rhs data with stimulus current
+        "intan/test_fcs_dc_250327_154333/info.rhs", # this is an example of only having dc amp rather than amp files
+        #"intan/test_fpc_stim_250327_151617/info.rhs", # wrong files Heberto will fix
     ]
 
     def test_annotations(self):

From 5256cd06b6e22232cfc8b66044cfbcf83bbe3d30 Mon Sep 17 00:00:00 2001
From: Heberto Mayorquin <h.mayorquin@gmail.com>
Date: Tue, 1 Apr 2025 13:29:53 -0600
Subject: [PATCH 19/20] clean up logic

---
 neo/rawio/intanrawio.py | 71 ++++++++++++++++++++++-------------------
 1 file changed, 39 insertions(+), 32 deletions(-)

diff --git a/neo/rawio/intanrawio.py b/neo/rawio/intanrawio.py
index 6bd466d43..3c393b2d8 100644
--- a/neo/rawio/intanrawio.py
+++ b/neo/rawio/intanrawio.py
@@ -899,19 +899,18 @@ def read_rhs(filename, file_format: str):
         names_to_count = [name for name in stream_names if name not in special_cases_for_counting]
         channel_number_dict = {name: len(stream_name_to_channel_info_list[name]) for name in names_to_count}
 
-        # Both DC Amplifier and Stim streams have the same number of channels as the amplifier stream
-        # if using the `header-attached` or `one-file-per-signal` formats
-        # the amplifier data is stored in the same place in the header so no matter what the DC amp
-        # uses the same info as the RHS amp.
+        # Each DC amplifier channel has a corresponding RHS2000 amplifier channel
         channel_number_dict["DC Amplifier channel"] = channel_number_dict["RHS2000 amplifier channel"]
-        if file_format != "one-file-per-channel":
-            channel_number_dict["Stim channel"] = channel_number_dict["RHS2000 amplifier channel"]
-            raw_file_paths_dict = create_one_file_per_signal_dict_rhs(dirname=filename.parent)
-        else:
+        
+        if file_format == "one-file-per-channel":
+            # There is a way to switch off amplifier and only keep the DC amplifier,
+            # so we need to count the number of files we find instead of relying on the header.
             raw_file_paths_dict = create_one_file_per_channel_dict_rhs(dirname=filename.parent)
             channel_number_dict["Stim channel"] = len(raw_file_paths_dict["Stim channel"])
-            # but the user can shut off the normal amplifier and only save dc amplifier
+            # Moreover, even if the amplifier channels are on the header their files are dropped
             channel_number_dict["RHS2000 amplifier channel"] = len(raw_file_paths_dict["RHS2000 amplifier channel"])
+        else:
+            channel_number_dict["Stim channel"] = channel_number_dict["RHS2000 amplifier channel"]
 
         header_size = f.tell()
 
@@ -925,7 +924,7 @@ def read_rhs(filename, file_format: str):
         memmap_data_dtype["timestamp"] = "int32"
         channel_number_dict["timestamp"] = 1
 
-    if file_format != "one-file-per-channel" or channel_number_dict["RHS2000 amplifier channel"] > 0:
+    if channel_number_dict["RHS2000 amplifier channel"] > 0:
         for chan_info in stream_name_to_channel_info_list["RHS2000 amplifier channel"]:
             chan_info["sampling_rate"] = sr
             chan_info["units"] = "uV"
@@ -967,29 +966,37 @@ def read_rhs(filename, file_format: str):
     # I can't seem to get stim files to generate for one-file-per-channel
     # so ideally at some point we need test data to confirm this is true
     # based on what Heberto and I read in the docs
+
+    # Add stim channels
     for chan_info in stream_name_to_channel_info_list["RHS2000 amplifier channel"]:
-        # we see which stim were activated
-        if file_format != "one-file-per-channel" or any(
-            [chan_info["native_channel_name"] in stim_file.stem for stim_file in raw_file_paths_dict["Stim channel"]]
-        ):
-            chan_info_stim = dict(chan_info)
-            name = chan_info["native_channel_name"]
-            chan_info_stim["native_channel_name"] = name + "_STIM"
-            chan_info_stim["sampling_rate"] = sr
-            # stim channel are complicated because they are coded
-            # with bits, they do not fit the gain/offset rawio strategy
-            chan_info_stim["units"] = "A"  # Amps
-            chan_info_stim["gain"] = global_info["stim_step_size"]
-            chan_info_stim["offset"] = 0.0
-            chan_info_stim["signal_type"] = 11  # put it in another group
-            chan_info_stim["dtype"] = "int16"  # this change is due to bit decoding see note below
-            ordered_channel_info.append(chan_info_stim)
-            # Note that the data on disk is uint16 but the data is
-            # then decoded as int16 so the chan_info is int16
-            if file_format == "header-attached":
-                memmap_data_dtype += [(name + "_STIM", "uint16", BLOCK_SIZE)]
-            else:
-                memmap_data_dtype["Stim channel"] = "uint16"
+        # stim channels are not always present in the header
+
+        if file_format == "one-file-per-channel":
+            # Some amplifier channels don't have a corresponding stim channel,
+            # so we need to make sure we don't add channel info for stim channels that don't exist.
+            # In this case, if the stim channel has no data, there won't be a file for it.
+            stim_file_paths = raw_file_paths_dict["Stim channel"]
+            amplifier_native_name = chan_info["native_channel_name"]
+            stim_file_exists = any([amplifier_native_name in stim_file.stem for stim_file in stim_file_paths])
+            if not stim_file_exists:
+                continue
+
+        chan_info_stim = dict(chan_info)
+        name = chan_info["native_channel_name"]
+        chan_info_stim["native_channel_name"] = name + "_STIM"
+        chan_info_stim["sampling_rate"] = sr
+        chan_info_stim["units"] = "A"  # Amps
+        chan_info_stim["gain"] = global_info["stim_step_size"]
+        chan_info_stim["offset"] = 0.0
+        chan_info_stim["signal_type"] = 11  # put it in another group
+        chan_info_stim["dtype"] = "int16"  # this change is due to bit decoding see note below
+        ordered_channel_info.append(chan_info_stim)
+        # Note that the data on disk is uint16 but the data is
+        # then decoded as int16 so the chan_info is int16
+        if file_format == "header-attached":
+            memmap_data_dtype += [(name + "_STIM", "uint16", BLOCK_SIZE)]
+        else:
+            memmap_data_dtype["Stim channel"] = "uint16"
 
     # No supply or aux for rhs files (ie no stream_id 1 and 2)
     # We have an error above that requests test files to help if the spec is changed

From 7128e6fe98d0e88499afed8ef17c1ac2ed5bb2ae Mon Sep 17 00:00:00 2001
From: Zach McKenzie <92116279+zm711@users.noreply.github.com>
Date: Wed, 2 Apr 2025 08:24:10 -0400
Subject: [PATCH 20/20] add more info to comments

---
 neo/rawio/intanrawio.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/neo/rawio/intanrawio.py b/neo/rawio/intanrawio.py
index ef6e89386..be665cb51 100644
--- a/neo/rawio/intanrawio.py
+++ b/neo/rawio/intanrawio.py
@@ -903,11 +903,11 @@ def read_rhs(filename, file_format: str):
         channel_number_dict["DC Amplifier channel"] = channel_number_dict["RHS2000 amplifier channel"]
         
         if file_format == "one-file-per-channel":
-            # There is a way to switch off amplifier and only keep the DC amplifier,
+            # There is a way to shut off saving amplifier data and only keeping the DC amplifier or shutting off all amplifier file saving,
             # so we need to count the number of files we find instead of relying on the header.
             raw_file_paths_dict = create_one_file_per_channel_dict_rhs(dirname=filename.parent)
             channel_number_dict["Stim channel"] = len(raw_file_paths_dict["Stim channel"])
-            # Moreover, even if the amplifier channels are on the header their files are dropped
+            # Moreover, even if the amplifier channels are in the header their files are dropped
             channel_number_dict["RHS2000 amplifier channel"] = len(raw_file_paths_dict["RHS2000 amplifier channel"])
         else:
             channel_number_dict["Stim channel"] = channel_number_dict["RHS2000 amplifier channel"]
@@ -970,7 +970,7 @@ def read_rhs(filename, file_format: str):
 
     # Add stim channels
     for chan_info in stream_name_to_channel_info_list["RHS2000 amplifier channel"]:
-        # stim channels are not always present in the header
+        # stim channel presence is not indicated in the header so for some formats each amplifier channel has a stim channel, but for other formats this isn't the case.
 
         if file_format == "one-file-per-channel":
             # Some amplifier channels don't have a corresponding stim channel,