load testfft 
windowtime=5;
samplerate=125;

%hardcode LFHP value
lfhp=0.75; %set the lo freq high pass to 0.75 hz.  Pulse rates below 45 bmp
%will not be recognized.  This is to help separate modulation due to heart
%pulse from modulation due to breathing rate

%hardecode maxhr (max heart rate) value
maxhr=3; %this is the max HR (in Hz) that will be used to determine the 
%fundamental heart pulse rate.  Fundamental pulse rates above 180bmp will
%not be recognized.

%define the number of spectra samples per window to go into fft
windowsamples=double(int32((windowtime*samplerate)));

%define the fft frequecies for the given fft size (window) and samplerate
freq=(0:floor(windowsamples/2)-1)'/windowsamples*samplerate;

%define the principal channel (channel with highest signal) to be used to
%determine the fundamental pulse rate and the HFLP value

principalchannel=find(mean(sd)==max(mean(sd)));

pcfft=fft(sd(1:windowsamples,principalchannel));
pcfft=pcfft(1:floor(windowsamples/2));
pcfft=pcfft.*conj(pcfft);

% find the highest freq ind that will be computed by the fft that is
% greater than or equal to the lfhp frequency (typically 0.75Hz)

startfreqind=find(freq>=lfhp,1,'first');

%find the lowest freq in that will be computed by the fft that is less than
%or equal to the max heart rate (typicall 180bpm or 3Hz).
maxhrfreqind=find(freq<=maxhr,1,'last');

%find the frequency that corresponds to the fft index with the highest
%amplitude between the lfhp and maxhr frequencies. This is the fundamental
%heart pulse frequncy for the principal channel

fundfreqind=find(pcfft==max(pcfft(startfreqind:maxhrfreqind)));

%define an arbitray "very high" frequency limit. this can be used to cut
%off anomlies in the data that have very hi frequencies or to set ifft
%components at very hi frequencies to zero real and imag. Above 8Hz will be
%considered very hi, fundamental and harmonics above this can be rejected
%(480 bmp)

hihifreq=find(freq<=8,1,'first');

%The signal modulation due to the heart puluse is not a pure sine or
%cosine, so the modulated signal has multiple harmonoics of the fundamental
%frequency.  The fundamental and first harmonic usually make up a +/-/+
%amplitude of the complex fft.  This range can be modified (ie extended to
%higher harmonics) if a more accurate modulated ifft is needed.

% hflp=2*freq(fundfreqind)+lfhp;
hflp=2*freq(fundfreqind);

stopfreqind=find(freq<=hflp,1,'last');

%Define the logical vector of fft indicies that are NOT the DC and in-band
%components

fftnot=true(windowsamples,1);
fftnot([1,startfreqind:stopfreqind,(end-stopfreqind+2):(end-startfreqind+2)])=0;
ffthinot=true(windowsamples,1);
ffthinot([1:stopfreqind,(end-stopfreqind+2):end])=0;
fftlonot=true(windowsamples,1);
fftlonot([1,startfreqind:end-startfreqind+2])=0;
ffthihinot=true(windowsamples,1);
ffthihinot([1:41,(end-41+2):end])=0;

bpFilt = designfilt('bandpassfir', 'FilterOrder', 25, 'CutoffFrequency1', .01, 'CutoffFrequency2', 8,'SampleRate', 25);
bpFilt2 = designfilt('bandpassfir', 'FilterOrder', 25, 'CutoffFrequency1', .01, 'CutoffFrequency2', 3,'SampleRate', 125);
bpFilt3 = designfilt('bandpassfir', 'FilterOrder', 25, 'CutoffFrequency1', .75, 'CutoffFrequency2', 2,'SampleRate', 125);

% f=fft(filtfilt(bpFilt,sd(1:625,principalchannel))/(sum(bpFilt.Coefficients)^2));
f=fft(filtfilt(bpFilt,sd(1:windowsamples,:))/(sum(bpFilt.Coefficients)^2));
ff=f;
f(fftnot,:)=0;%all non interesting frequencies real and imag set to zero
ff(fftnot,:)=0+1i*imag(ff(fftnot,:));%all non interesting frequencies real only set to zero
fff=ff;
fff(ffthihinot,:)=0;%same as ff, but freqs above 8hz imag set to zero

% plot(linspace(0,5,625),[sd(1:625,principalchannel),ifft(f),ifft(ff),ifft(fff)]);grid
% legend('Signal','All Zero','Imag Kept','Imag Only Below 8Hz');shg