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)));

%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=logical(ones(windowsamples,1));
fftnot([1,[startfreqind:stopfreqind],[(end-stopfreqind+2):(end-startfreqind+2)]])=0;
ffthinot=logical(ones(windowsamples,1));
ffthinot([[1:stopfreqind],[(end-stopfreqind+2):end]])=0;
fftlonot=logical(ones(windowsamples,1));
fftlonot([1,[startfreqind:end-startfreqind+2]])=0;
ffthihinot=logical(ones(windowsamples,1));
ffthihinot([[1:41],[(end-41+2):end]])=0;

f=fft(sd(1:625,principalchannel));
ff=f;
f(fftnot)=0;%all non interesting frequencies real and imag set to zero
ff(fftnot)=0+1*i*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,20),ifft(f),ifft(ff),ifft(fff)]);grid
legend('Signal','All Zero','Imag Kept','Imag Only Below 8Hz');shg