function [pdsabs,pdsabs2,sdfl,sdflh]=signaltopds(sd,samplerate,windowtime,filterorder,loband,hiband,filteron)

%USAGE:
%[pdsabs,sdfl,sdflh]=signaltopds(sd,samplerate,windowtime,filterorder,loband,hiband)
%
%Bandpass filters the spectral data set and returns the pulse differential
%spectra along with the lopass filtered data and the bandpass filtered data
%
%sd is the dark corrected spectra
%sample rate is the rate the spectra are collected
%window time is the width of the window (in seconds) to move over sd and
%   perform ffts
%filterorder if the order of the for the lo and hi pass butterworth
%loband and hiband are the lower and upper bandpass edges (in Hz) for the
%   butterworth
%
%pdsabs is the pulse differential spectra
%sdfl is the lopass filtered spectra
%sdflh is the bandpass filtered spectra
%
%For example, if 30 seconds of spectra are collected with 1900μs of
%integration time and 600μs of readout (for a total of 2500μs spectral
%rate), there would be 12000 spectra. Using a window of 10.24 seconds and a
% bandpass betwee 0.5Hz and 7.5Hz:
%
%[pdsabs,sdfl,sdflh]=signaltopds(sd,400,10.24,6,0.5,7.5)
%
% loband=0.5;
% hiband=7.5;
% filterorder=6;
% samplerate=400;
% windowtime=10.24;
%
%USE [pdsabs,~,~]=signaltopds(sd,400,10.24,6,0.5,7.5)
%To return only the pds spectra

%define a the size of the window in samples based on the with of the window
%in seconds
windowsamples=(windowtime*samplerate);
%define the fft frequecies for the given fft size (window) and samplerate
freq=[0:windowsamples/2-1]'/windowsamples*samplerate;
%find the highest frequency that will be computed by the fft that is less
%than or equal to the loband frequency
startfreq=find(freq<=loband,1,'last');
stopfreq=find(freq<=hiband,1,'last');
%compute the normalized freqeuncy for the lo pass filter. It should only
%pass frequencies lower than the hiband cutoff frequency. Note that
%highband cannot exceed the sample rate/(4pi) since the butterworth
%normalized filter frequency must be between 0 and 1
flo=2*pi*hiband/(samplerate/2);
%compute the normalized freqeuncy for the hi pass filter. It should only
%pass frequencies higher than the loband cutoff frequency.
fhi=2*pi*loband/(samplerate/2);
%compute the butterworth filter coeffecients
[a,b]=butter(filterorder,flo,'low');
[aa,bb]=butter(filterorder,fhi,'high');
%filter the signal with the lo pass filter
sdfl=filtfilt(a,b,sd);
%filter the signal with the hi pass filter
sdflh=filtfilt(aa,bb,sdfl);

%run a sliding window of "windosamples" across the sample set. In each
%window compute the fft on each channle in the window. Use interpolation to
%find the amplitude of the fundamental frequency, 1st harmonic and 2nd
%harmonics for each channel.  Assume a peak width +/- 1/3 the value of the
%fundamental frequency. This should be sufficiently far away from the
%loband frequency and allow sufficient interpolation across the peak to 
%find a peak even when the peak is between fft frequencies

%pre allocate the pdsabs variable so it doesn't slow down the loop
pdsabs=zeros(size(sd,1)-windowsamples+1,size(sd,2));
pdsabs2=pdsabs;

parfor i=1:size(sd,1)-windowsamples+1
% for i=1:size(sd,1)-windowsamples+1
    %compute the fft on the windos of the hi and lo pass filtered signal
    if filteron
        fftamp=abs(fft(sdflh(i:i+(windowsamples-1),:))/windowsamples);
    else
        fftamp=abs(fft(sd(i:i+(windowsamples-1),:))/windowsamples);
    end

    %It is assumed that pixel 36 (69 for trans tune) has the strongest reflectance signal the 
    %data set making the pulse signal to noise highest.  Using this pixle
    %to find the fft element the fundmental pulse frequency occurs.  Do the
    %same for the 1st and 2nd harmonic of the pulse frequency.
    fundpix=69;
    fundid=find(fftamp(1:floor(windowsamples/2),fundpix)==max(fftamp(startfreq:stopfreq,fundpix)));
    %the first fft element is the DC (mean) value.
    freqid=fundid-1;
    %define the 1/3 width value of the fundamental pulse frequency
    sig=floor((freqid)/3);    
    firstid=find(fftamp(1:floor(windowsamples/2),fundpix)==max(fftamp(freqid*2+1-sig:freqid*2+1+sig,fundpix)));
    secondid=find(fftamp(1:floor(windowsamples/2),fundpix)==max(fftamp(freqid*3+1-sig:freqid*3+1+sig,fundpix)));
    %for the fft element where the fundamental, 1st and 2nd harmonic
    %frequencies occur, interpolate the fft with 100 points between
    %elements.  Ideally, this would be a weighted spline fit, but for
    %computational economy, spline or cubic are sufficient
    out=interpn(fundid-sig:secondid+sig,fftamp(fundid-sig:secondid+sig,1:256),fundid-sig:.01:secondid+sig,'cubic')';
    %setup a 2 col matrix of the 100x resolution pixel 36 data to find the
    %elements for the fundamental, 1st and 2nd harmonics for the hi
    %resoltuion interpolated fit
    peaksig=[[fundid-sig:.01:secondid+sig]',out(:,fundpix)];
    %find the elements of the hi resoltuion fit data
    fundidhr=find(peaksig(:,2)==max(peaksig(find(peaksig(:,1)==fundid)-100*sig:find(peaksig(:,1)==fundid)+100*sig,2)));
    firstidhr=find(peaksig(:,2)==max(peaksig(find(peaksig(:,1)==firstid)-100*sig:find(peaksig(:,1)==firstid)+100*sig,2)));
    secondidhr=find(peaksig(:,2)==max(peaksig(find(peaksig(:,1)==secondid)-100*sig:find(peaksig(:,1)==secondid)+100*sig,2)));
    hrids=[fundidhr firstidhr secondidhr]';
    %define the PDS absorbance using the fundamental, 1st and 2nd harmonic
    %frequency amplitudes and the DC value of the window.  The PDS spectrum
    %is defined as PDS = log10(1 + ΔI/I), where I is the DC intensitiy
    %value of the channel in the window and ΔI is the peak to peak
    %intensity of the channel due to pulse.  ΔI is approximated by the sum
    %of the amplitudes for the fundamental, 1st and 2nd harmonic
    %frequencies of the pulse.  Summing more harmoincs is recommended if
    %the SNR allows.  The FFT amplitude is multiplied by 2 since the
    %amplitude of each frequency is split between the computed frequencies
    %and the mirroed aliases.  Since the FFT amplitude represents the
    %magnitude of the deviation of the signal from the DC value an
    %additional 2x is required to match the peak to peak of ΔI.
    %COMMENT: PDS assumes the lowest signal of the pulse train is the DC
    %value and not the mean of the pulse. It may be necessary to subtract
    %1/d of ΔI from the DC value when calculating the PDS spectrum as shown
    %in the next line.  NOTE, computing the mean in this line is needed if
    %the filtered data is used since the mean is removed.
    %NOTE: in the 1450nm - 1550nm range, the transmissive signal may be
    %close to the noise floor.  If the filtered data does not produce an
    %FFT amplitude that is realistic (check the fft in this range to see if
    %the amplitudes are above the fft noise floor), it may be necessary to
    %turn off the filtering and extend the number of samples used in the
    %fft.  For example, a 30s data set sampled at 400hz (hi gain,
    %transmission) may not have significant fft values in this low signal 
    % range.  While the rest of the range may have high snr pds amplitudes
    % with only 5 or 10 s of data (2000 or 4000 spectra), turning off
    % filtering and extending the number of spectra to 9000 or 10000 (22.5
    % to 25s of data) may increase the snr of the pds amplitudes in the low
    % signal range.  Caution should be used when extending the number of
    % spectra (increasing the sample time input the fft) as cardiac pulse
    % rate may be less consistent with longer sample duration.

    pdsabs(i,:)=log10(1+4*sum(out(hrids,:))./(mean(sd(i:i+(windowsamples-1),:))-2*sum(out(hrids,:))));
    pdsabs2(i,:)=log10(1+4*sum(out(hrids,:))./mean(sd(i:i+(windowsamples-1),:)));
end
