function [time, B_output, x, f, t, spectraltime_axis, fs, titlestr, siteData, mlt_time, mlt_starttime, xspectrstart, spectraltime_axis_mlt] = carisma_icm_psd(folder, starthour, stophour, npt, direction, filterchoice, utc_or_mlt)

global overlap
mlt_time = [];
mlt_starttime = [];
xspectrstart = [];

[time,B_Hunfilt,B_Dunfilt,Qual,siteData] = importCARISMAICMFolder(folder); % Import ICM data

time.TimeZone = 'UTC';

fs = siteData.Fs; % sampling frequency, usually 20 or 100

if strcmp(filterchoice,'yes')
%     [pH,sH,muH] = polyfit((1:numel(B_Hunfilt))',B_Hunfilt,20);
%     f_H = polyval(pH,(1:numel(B_Hunfilt))',[],muH);
%     B_H = B_Hunfilt - f_H;
%     [pD,sD,muD] = polyfit((1:numel(B_Dunfilt))',B_Dunfilt,6);
%     f_D = polyval(pD,(1:numel(B_Dunfilt))',[],muD);
%     B_D = B_Dunfilt - f_D;
    if fs == 100
        wn_high=0.3/(fs/2); % high-pass filter.  Cannot remember why 2.8*0.3 = 0.84.  IAR is 0.5-3Hz.  Trying 0.5 (July 16,2018) instead. It breaks.  It does not filter out all the DC, so the signal doesn't sit at zero and that breaks my code, because it needs everything to be flat.
        wn_low = 15/(fs/2);
    else
        wn_high=0.3/(fs/2); % high-pass filter
        wn_low = 0.9;
    end
    filter_coeff_high=fir1(fs*20,wn_high,'high',hanning(fs*20+1)); % filter coefficient
    filter_coeff_band = fir1(fs*20,[wn_high wn_low],'',hanning(fs*20+1)); % filter coefficient
    B_H = filtfilt(filter_coeff_high,1,B_Hunfilt); % filter the data
    B_D = filtfilt(filter_coeff_high,1,B_Dunfilt);
    B_H_band = filtfilt(filter_coeff_band,1,B_Hunfilt); % filter the data
    B_D_band = filtfilt(filter_coeff_band,1,B_Dunfilt);
end

% Create time vector bounds in seconds
abs_start = datetime(year(time(1)), month(time(1)), day(time(1)),starthour,0,0);
abs_start.TimeZone = 'UTC';
abs_stop = datetime(year(time(end)), month(time(end)), day(time(end)),stophour,59,59.95);
abs_stop.TimeZone = 'UTC';

% Truncate data
B_H=B_H(time >= abs_start & time <= abs_stop);
B_D=B_D(time >= abs_start & time <= abs_stop);
B_H_band=B_H_band(time >= abs_start & time <= abs_stop);
B_D_band=B_D_band(time >= abs_start & time <= abs_stop);
Qual = Qual(time >= abs_start & time <= abs_stop);
time = time(time >= abs_start & time <= abs_stop);

% Generate titlestr
titlestr = strcat('CARISMA ICM',{' '},direction,'-component:',{' '}, datestr(time(1),'yyyy-mm-dd'), {' '}, siteData.siteID);

%----------------------
% specify sliding window fft parameters
overlap=round(0.3*npt);

% Calculate psd parameters, x (spectrogram), f (frequency), t (time)
if direction == 'H'
    [x, f, t]=spectrogram(B_H, hanning(npt), overlap, npt, fs); % fft of B_H
    B_output = B_H;
elseif direction == 'D'
    [x, f, t]=spectrogram(B_D, hanning(npt), overlap, npt, fs); % fft of B_D
    B_output = B_D;
elseif direction == 'H\_band'
    [x, f, t]=spectrogram(B_H_band, hanning(npt), overlap, npt, fs); % fft of B_H
    B_output = B_H_band;
elseif direction == 'D\_band'
    [x, f, t]=spectrogram(B_H_band, hanning(npt), overlap, npt, fs); % fft of B_H
    B_output = B_H_band;
end

if strcmp(utc_or_mlt,'MLT') % FIXMEEEEEEEEEE
    switch siteData.siteID
        case 'DAWS'
            mlt_time = datetime(time,'TimeZone','America/Los_Angeles','Format','dd-MMM-uuuu HH:mm:ss.SS Z');
        case {'MSTK','FSMI'}
            mlt_time = datetime(time,'TimeZone','America/Denver','Format','dd-MMM-uuuu HH:mm:ss.SS Z');
        case {'FCHU','PINA','RABB','THRF'}
            mlt_time = datetime(time,'TimeZone','America/Chicago','Format','dd-MMM-uuuu HH:mm:ss.SS Z');     
    end
    utc_time = time; % save time as utc_time
    mlt_starttime = find(datenum(mlt_time) == datenum(utc_time(1))); % find lower truncation limit
    %time = mlt_time(mlt_starttime:end); % truncate mlt_time and set it as "time"
    %B_output = B_output(mlt_starttime:end); % truncate B_output
    Qual = Qual(mlt_starttime:end); % truncate Qual
    
    % truncate x
    trunc_ratio = mlt_starttime/length(mlt_time); % find ratio of truncated data in mlt_time
    szx = size(x); 
    xspectrstart = round(trunc_ratio*szx(2)); % define start point for truncation
    %x = x(:,spectrstart:end); % truncate!
    abs_start = time(1); % identify upper truncation limit for spectral data
    abs_stop = time(end); % identify upper truncation limit for spectral data
    
end

% generate time-axis for time series
%sz1=size(time); rawtime_axis=linspace(hour(time(1)), hour(time(end))+1, sz1(1));

% generate time-axis for spectral plot
sz2=size(x); N = 1:sz2(2); midwinpoint = npt/2+(N-1)*(npt-overlap); 
spectraltime_axis = time(midwinpoint);%spectraltime_axis=linspace(hour(abs_start), hour(abs_stop)+1, sz2(2));
if strcmp(utc_or_mlt,'MLT')
spectraltime_axis_mlt = mlt_time(midwinpoint);
end
end
