function PSP_finder(fina,pspwavename,length_PSP,pspthresh); 

%fina :  the long recording composed of many PSPs/PSCs to be analyzed

%pscwavename : the representative risetime waveform to be used in finding PSPs/PSCs 

%length_PSP :  for example, if length_PSP = 10, then this 
%   program would assume that the wave to be analyzed is composed of PSCs or
%   PSPs that are approximately 10 msec in duration.

%pspthresh :  a somewhat arbitrary threshold set for determining whether a
%   correlation match is significant or not.  A trial and error approach might
%   be warranted here.  Usually pspthresh is greater than 90

% this program assumes both pscwavename and fina are matfiles

oldwidth=0;

pscwavename=pspwavename;

tfirst=0; %ycount=0; yyymem=[];
pscwave=load(pscwavename);     
pscwave=getfield(pscwave,pscwavename);  

increment = length (pscwave);  %pscwave is an estimate of the average risetime for the PSCs/PSPs

tic;
ttttc=[];
ttt=1;
aaaa=0;
cccc=0;


increment=round(length_PSP/5);   %increment = 3;
qcount=0;detour=1;
suffix=num2str(pspthresh);
pspthresh=pspthresh/100;
savename = [fina,'_psptimes',suffix];
    
%----------------  LowPass filter --------------

namefile=fina;                          %if you are using .mat files instead of
axowave=load(fina); 
axowave=getfield(axowave,fina);
axowave_old=axowave;
mean_axo=mean(axowave);
axowave=axowave-mean(axowave);
normalized(1:100000)=0;
lenaxo=size(axowave,1);
normalized(100001:(size(axowave,1)+100000))=axowave;
normalized((size(axowave,1)+100000):(size(axowave,1)+200000))=0;
Sample_rate_Hz = 1000;
lowpass_Hz = 1;
Wn=(lowpass_Hz/(Sample_rate_Hz/2));
b = fir1(1000,Wn,'low');
%freqz(b,1,500,1000);
threshwave = filter(b,1,normalized);
normalized=normalized';
threshwave=threshwave';
cutsize=Sample_rate_Hz/2;
normalized((lenaxo+100000):(size(normalized,1)))=[];
normalized(1:100000)=[];
threshwave((lenaxo+100000+cutsize):(size(threshwave,1)))=[];
threshwave(1:(100000+(cutsize-1)))=[];
threshwave=threshwave+mean_axo;
normalized=normalized+mean_axo;

%threshwave=threshwave-3;
threshwave_old = threshwave;


milepost='lowpass filter'

%---------------  findPSCs --------------------------

%function [tempwins]=findpscs_1khz_func(fina,finb,fint,increment) 
SampleRate = 1000;     % points/second
%increment = 16;     % [points] - template length
threshold = 0.1;
starttime = 1;
resolutionvar=increment;
width = increment/resolutionvar;    % [points] - step size
wincount = 0;
counter=1;



pscnpnts = length(pscwave);
axonpnts = length(axowave);


num = max(pscnpnts,axonpnts)-increment;



winnerlist = [];
tempwins = [];

endchunk = ceil((axonpnts-(2*SampleRate))/SampleRate);
winners_counter = 1;  counter=1;
        delta_psc = pscwave;
        
for chunk=1:endchunk

%dagwave=load(finb);     
%dagwave=getfield(dagwave,finb);
axowave = axowave_old;
axowave(((chunk+1)*SampleRate+1):axonpnts)=[];
%threshwave=load(fint);
%threshwave=getfield(threshwave,fint);
threshwave=threshwave_old;
threshwave(((chunk+1)*SampleRate+1):axonpnts)=[];
%counter
if (chunk>1)
    axowave(1:(chunk*SampleRate))=[];
    threshwave(1:(chunk*SampleRate))=[];
end;
    % while (counter < num)
    while (length(axowave)>pscnpnts)
        %delta_axo = axowave(counter:counter+increment);
        delta_axo = axowave(1:increment);

         % using > or < depends on vclamp vs. current clamp in next line   
         if ((min(delta_axo)>threshwave(1)) && (delta_axo(increment)>delta_axo(1)))
               
              
               shift_CCG =  xcov(delta_psc,delta_axo,'coeff');
               
               if (shift_CCG(increment)>pspthresh) 
                    tempwins(winners_counter,1) = shift_CCG(increment);
                    tempwins(winners_counter,2) = counter;
                    winners_counter = winners_counter+1;
               end;    

         end;
         counter=counter+1;
         axowave(1)=[];
       %  threshwave(1)=[];


    %       SDpsc = 0;
    %       SDaxo = 0;
    %       shift_CCG = [];


  end;

      counter=counter+increment;


end;

%namefile=strrep(namefile,'.txt','');
%result = [namefile,'psctimes'];


%save (result,'tempwins','-ASCII');

pscsfound = tempwins;

%clear




milepost='PSPfinder'

%--------------------------- end PSCfinder-------------------------


%----------------------------Begin IdIPSPs-------------------------

%Macro  idpsps(pspwave,correlwave,timewave,threshname,thresh)
%String pspwave; String correlwave; String timewave; String threshname; Variable thresh
%Silent 1;PauseUpdate
axowave=axowave_old;
threshwave = threshwave_old;
thresh = pspthresh;
correltrim = pscsfound(:,1);
timewave_uncut = pscsfound(:,2);
timewave = timewave_uncut; 


wavelength=ceil(length_PSP);    %wavelength=ceil(length_PSP/2);

xxx=1; yyy=1;
while (length(timewave)>xxx)

	if ((timewave(xxx+1)-timewave(xxx)) < wavelength)
		if (correltrim(xxx+1) > correltrim(xxx))
			correltrim(xxx)=[];
			timewave(xxx)=[];
			if (xxx>0)
                xxx=xxx-1;
            end;

			else			
				correltrim(xxx+1)=[];
				timewave(xxx+1)=[];
				if (xxx>0)
                    xxx=xxx-1;
                end;

       end;
    end;
    
    if (xxx>0)
        if (correltrim(xxx) < thresh)
        	correltrim(xxx)=[];
            timewave(xxx)=[];
            if (xxx>0)
                xxx=xxx-1;
            end;

        end;
    end;
	%WaveStats/Q/R=((timewave[xxx]),(timewave[xxx]+30)) $pspwave
    if (xxx>0)
        V_max = max(axowave(timewave(xxx):timewave(xxx)+wavelength));
        if (V_max<threshwave(timewave(xxx)))
            correltrim(xxx)=[];
            timewave(xxx)=[];
            if (xxx>0)
            	xxx=xxx-1;
            end;
            yespass = 0;
        end;
    end;

    	xxx=xxx+1;

end;    
%duplicate/o correltrim tencorrel
%tencorrel*=10


%------------------- end IDpsps --------------------

milepost='IDpsps'

%length(timewave)

save(savename,'timewave');






    



%save (namefile, 'tempwins','-ASCII');  %can be turned off.  leaving it on can be
                                     %useful for diagnostics

 



