FPGA

using Pkg
Pkg.activate("/home/pi/.julia/environments/v1.9")  # Replace with your Julia environment path

using RedPitayaDAQServer
using CSV
using DataFrames
using DSP          # For signal processing functions
using Statistics   # For statistical functions like mean, std
using Dates        # For timestamping
# using Plots       # Uncomment if you wish to add visualization (feature 14)

# Set the Red Pitaya IP and port 5025
rp = RedPitaya("rp-f0ac8a.local", 5025)

# Set server in CONFIGURATION mode to prepare signal acquisition
serverMode!(rp, CONFIGURATION)

# Set decimation to 1250 for a sampling rate of 100 kS/s
decimation!(rp, 1250)  # 125 MS/s / 1250 = 100 kS/s

# Set number of samples per period and frame
samples_per_period = 16384  # Adjust if needed
periods_per_frame = 1       # For continuous acquisition

samplesPerPeriod!(rp, samples_per_period)
periodsPerFrame!(rp, periods_per_frame)
triggerMode!(rp, INTERNAL)

# Start acquisition
serverMode!(rp, ACQUISITION)
masterTrigger!(rp, true)

# Initialize counters
global samples_received = 0
global frames_received = 0
global neutron_events_count = 0
global neutron_events_in_last_second = 0
global last_cps_update_time = 0.0

# Initialize background level for adaptive thresholding (feature 10) and background subtraction (feature 9)
global background_level = 0.0
global noise_std = 0.0

# Function to detect and discriminate neutron events
function isNeutronEvent(signal, sampling_rate)
    # Ensure the signal is a 1D array
    signal = vec(signal)

    # Signal Smoothing and Noise Reduction (Feature 4)
    # Noise filtering using a low-pass filter
    cutoff_freq = 40_000      # 40 kHz cutoff frequency

    # Normalize cutoff frequency
    nyquist_freq = sampling_rate / 2
    normalized_cutoff = cutoff_freq / nyquist_freq

    # Ensure normalized_cutoff is less than 1
    if normalized_cutoff >= 1.0
        error("Normalized cutoff frequency must be less than 1")
    end

    # Design a Butterworth low-pass filter
    filter_order = 5
    butter_filter = digitalfilter(Lowpass(normalized_cutoff), Butterworth(filter_order))

    # Apply the filter
    filtered_signal = filtfilt(butter_filter, signal)

    # Baseline correction
    baseline = mean(filtered_signal[1:10000])  # Use first 10000 samples
    adjusted_signal = filtered_signal .- baseline

    # Background Radiation Subtraction (Feature 9)
    adjusted_signal .= adjusted_signal .- background_level

    # Adaptive Thresholding (Feature 10)
    # Update noise standard deviation
    global noise_std = std(adjusted_signal[1:10000])  # Use first 10000 samples
    pulse_threshold = background_level + 5 * noise_std  # Threshold set to 5 sigma above background

    # Detect pulses using custom peak detection
    # Identify indices where the signal crosses the threshold from below
    above_threshold = adjusted_signal .> pulse_threshold
    crossing_indices = findall(diff(above_threshold) .== 1) .+ 1  # Indices where signal rises above threshold

    # Initialize counter
    neutron_events_detected = 0

    # Analyze pulses
    for idx in crossing_indices
        # Define a window around the crossing point for pulse analysis
        window_size = 200  # Adjust based on expected pulse width
        start_idx = idx
        end_idx = min(idx + window_size - 1, length(adjusted_signal))
        pulse = adjusted_signal[start_idx:end_idx]

        # Check if pulse peak is above threshold
        pulse_amplitude = maximum(pulse)
        if pulse_amplitude < pulse_threshold
            continue  # Skip pulses that do not have a peak above the threshold
        end

        # Pulse Shape Discrimination (PSD) (Feature 1)
        fast_window_samples = min(50, length(pulse))  # Number of samples for fast component
        fast_component = sum(pulse[1:fast_window_samples])
        total_component = sum(pulse)
        psd_ratio = fast_component / total_component

        # Set PSD ratio threshold based on calibration
        psd_threshold = 0.3  # Adjust based on calibration

        # Decision logic based on PSD
        if psd_ratio < psd_threshold
            neutron_events_detected += 1

            # Event Logging and Time Stamping (Feature 12)
            event_time = now()
            # Log event data
            open("events_log.csv", "a") do io
                CSV.write(io, DataFrame(timestamp=Dates.format(event_time, "HH:MM:SS.sss"),
                                        psd_ratio=psd_ratio),
                          writeheader=false)
            end

            # For Visualization (Feature 14)
            # Uncomment and adjust the following lines if you wish to store pulse data for plotting
            # open("pulse_data.csv", "a") do io
            #     CSV.write(io, DataFrame(time=(start_idx:end_idx)/sampling_rate, amplitude=pulse), writeheader=false)
            # end
            # Continue processing other pulses
        end
    end

    return neutron_events_detected
end

# Track time
start_time = time()
global last_cps_update_time = start_time

# Create CSV file for CPS data
csv_file = "/home/pi/Desktop/Desktop/Neutron/rp_test/test.csv"
df = DataFrame(timestamp=Float64[], cps=Float64[])
CSV.write(csv_file, df)

# Sampling rate
sampling_rate = 125_000_000 / 1250  # 100 kS/s

# Main loop
while true
    fr = currentFrame(rp)
    uCurrentPeriod = readFrames(rp, fr, 1)

    # Ensure the signal is in the correct format
    signal = uCurrentPeriod[:]

    # Update counters
    global samples_received += samples_per_period * periods_per_frame
    global frames_received += 1

    # Update background level for background subtraction (Feature 9)
    # Assuming neutron events are rare, we can update background level continuously
    global background_level = mean(signal[1:10000])  # Use first 10000 samples

    # Check for neutron events
    neutron_events_detected = isNeutronEvent(signal, sampling_rate)

    # Update counters
    global neutron_events_count += neutron_events_detected
    global neutron_events_in_last_second += neutron_events_detected

    # Time calculations
    elapsed_time = time() - start_time
    samples_per_second = samples_received / elapsed_time
    frames_per_second = frames_received / elapsed_time

    # Update CPS every second
    current_time = time()
    if current_time - last_cps_update_time >= 1.0
        cps = neutron_events_in_last_second / (current_time - last_cps_update_time)

        # Display data
        println("Samples per second: $samples_per_second")
        println("Frames per second: $frames_per_second")
        println("Neutron Counts Per Second (CPS): $cps")

        # Write to CSV with timestamp (Feature 12)
        open(csv_file, "a") do io
            CSV.write(io, DataFrame(timestamp=current_time - start_time, cps=cps), writeheader=false)
        end

        # Reset counter
        global neutron_events_in_last_second = 0
        global last_cps_update_time = current_time
    end

    # Optional sleep to prevent high CPU usage
    # sleep(0.01)
end

# Stop acquisition (this part will not be reached in the infinite loop)
masterTrigger!(rp, false)
serverMode!(rp, CONFIGURATION)