function FluoData = DemodFluo(varargin)
    %     Function to create a read mat file output of Spike2 and generate
    %     Fluorescence data
    % 
    %     Input:   
    % 
    %     Output:
    %         FluoData: fluorescence data consist of following value
    %    Grady KaraZhang, March 2023
       
        %%
        p = inputParser;
        addParameter(p,'mod_freq',311, @isnumeric);
        addParameter(p,'Vscale',1/(2.0^15)*10, @isnumeric);
        addParameter(p,'fs',10000, @isnumeric);
        addParameter(p,'basepath',pwd, @isstr);
        addParameter(p,'basename',[], @isstr);
        addParameter(p,'Nfilter_baseline',3, @isnumeric);
        addParameter(p,'N_butter',4, @isnumeric);
        addParameter(p,'fcutoff',20, @isnumeric);
        parse(p,varargin{:});

        mod_freq = p.Results.mod_freq;
        Vscale = p.Results.Vscale;
        fs = p.Results.fs;
        basepath = p.Results.basepath;
        Nfilter_baseline = p.Results.Nfilter_baseline;
        N_butter = p.Results.N_butter;
        fcutoff = p.Results.fcutoff;
        if isempty(p.Results.basename)
            basename = bz_BasenameFromBasepath(basepath);
        else
            basename = p.Results.basename;
        end
     
        load([basename,'.mat'],'OXT_UF','OXT_Em','lfp_time','pulses');
    
        %
        %% Read unfiltered Oxt and Ach signal together with their modulation channels
        if ~isfield(OXT_UF,'ref')
            error('No reference signal found in the mat file')
        end
        if isa(OXT_UF.ref, 'double')
            disp('Fluo_ref is of type double.');
            Fluo_ref = OXT_UF.ref;
            Fluo_emission_unfiltered = OXT_UF.data;
        else
            disp('Fluo_ref is not of type double.');
            Fluo_ref = double(OXT_UF.ref)*Vscale;
            Fluo_emission_unfiltered = double(OXT_UF.data)*Vscale;
        end
        Vraw = double(OXT_Em.data)*Vscale;
        tt = OXT_Em.time;
        pulses_peak_times = mean(pulses.timestamps,2)-lfp_time(1);
        samplesToRemove = 1000;
        tt1 = downsample(tt,50);
        % Preallocate the peakLocation array
        peakLocation = zeros(size(pulses_peak_times));
        for i = 1:length(pulses_peak_times)
            [~, index] = min(abs(tt1 - pulses_peak_times(i)));
            peakLocation(i) = index*50;
                 % Calculate the range of samples to remove
            rangeStart = max(1, peakLocation(i) - samplesToRemove);
            rangeEnd = min(length(Vraw), peakLocation(i) + samplesToRemove);
            
            % Remove the samples within the defined range
            %Vraw(rangeStart:rangeEnd) = NaN;
        end
        %Vraw = fillmissing(Vraw, 'nearest');
        %%
        [b,a] = butter(N_butter,fcutoff/(fs/2));
        %% demodulate signals
        % if fmods are unknown can be found by checking the peak location of power spectra
        % will be used for band pass filtering but can be ignored also. To be
        % consistent with some other group I added. During demodulation we
        % will remove out of band signal.
        fmod_oxt = mod_freq;
        %probably extra filtering is not helpful here, can be removed
        Fluo_demod_emission_ = filtfilt(b,a,quadrature_demod(Fluo_emission_unfiltered, Fluo_ref,fs,fcutoff,fmod_oxt));
        Ndown = round(fs/100);
        Fluo_demod_emission = Fluo_demod_emission_(2*fs+1:end-2*fs);
        Fluo_demod_emission_ds = downsample(Fluo_demod_emission,Ndown);
        Ns = length(Fluo_demod_emission_ds);
        taxis = linspace(0,Ns/100,Ns);
        %% Calculate baseline
        [Fluo_dFF,Fluo_baseline] = calc_dFF(taxis,Fluo_demod_emission_ds,Nfilter_baseline,0);
        
        %% Generate output data
        FluoData.time = taxis';
        FluoData.Fluo_dFF = Fluo_dFF;
        FluoData.Fluo_baseline = Fluo_baseline;
        FluoData.Fluo_ds = Fluo_demod_emission_ds;
    
    end