Tue, 09/30/2008 - 03:05 pm
I am having a devil of a time coming up with a way to properly identify peak edges in my data. I'd like to use FindPeak, but I don't understand how it works very well. My data consists of a long wave with amplitude data- there are a series of "peaks" (actually 1-D MRI images), but they don't have any conventional shapes like a Gaussian or Lorentzian. I can get FindPeak to get the peak value right, but I have little to no interest in that. It's the edges that I'm after, and up till now, I've been identifying them by sight. This is no longer feasible. FindPeak has the nasty habit of telling me peak edges that are no where near correct. I think it is overly sensitive to what defines a peak. My data is a tad noisy, but what FindPeak is giving me is absolutely ridiculous. I am looking at a "peak" right now that is about 500 points wide, and Igor is telling me it's about 5 points wide. That makes no sense. It appears to be identifying amplitude swings from noise as actual peak edges. I don't understand the box size functionality, and just trying random numbers for the /b=value does not seem to help.
Does anyone have any suggestions on how to find peak edges? I am including an example of my data to better show what I am talking about (as both an image and an Igor binary). I'd visually identify points 2205 (roughly) and 2747 as the edges of this peak. This is a snippet of an actual dataset, which consists of a series of these guys, one right after the other. I'm sorry to ask such a basic question, but I've already invested a day on this, and am getting nowhere. I appreciate your help and feedback very much!
Peak1.ibw
Hope this helps! (Be sure to display w1 and w2 along with your original data to see what is going on.)
wave w; variable nsmooth
Duplicate /o w w1
Smooth nsmooth, w1
Differentiate w1 /D=w2
wavestats /q w2
print abs(V_minloc-V_maxloc)
// Killwaves w1,w2
End
John Bechhoefer
Department of Physics
Simon Fraser University
Burnaby, BC, Canada
September 30, 2008 at 06:56 pm - Permalink
V_npnts= 4096; V_numNaNs= 0; V_numINFs= 0; V_avg= 677.141;
V_Sum= 2.77357e+06; V_sdev= 1922.9; V_rms= 2038.42;
V_adev= 1218.66; V_skew= 2.5337; V_kurt= 4.69783;
V_minloc= 2139; V_maxloc= 2486; V_min= -118.567; V_max= 7024.83;
V_startRow= 0; V_endRow= 4095;
pulsestats/b=55 /f=0.99 /L=(V_max, peak1[V_endrow]) /M=(pnt2x(peak1, V_maxLoc) - pnt2x(peak1, V_minLoc)) Peak1
V_Flag= 1; V_PulseLoc1= 2218.27; V_PulseLoc2= 2750.07;
V_PulseLoc3= NaN; V_PulsePolarity= 1; V_PulseLvl0= 7024.83;
V_PulseLvl123= 681.587; V_PulseLvl4= -23.2181;
V_PulseAmp4_0= -7048.05; V_PulseWidth2_1= 531.797;
V_PulseWidth3_2= NaN; V_PulseWidth3_1= NaN;
V_Flag= 1; V_PulseLoc1= 2184.2; V_PulseLoc2= 2768.5;
V_PulseLoc3= NaN; V_PulsePolarity= 1; V_PulseLvl0= 7024.83;
V_PulseLvl123= 47.2623; V_PulseLvl4= -23.2181;
V_PulseAmp4_0= -7048.05; V_PulseWidth2_1= 584.295;
V_PulseWidth3_2= NaN; V_PulseWidth3_1= NaN;
Depending on what values you use for the /F flag, the V_PulseWidth2_1 value gets pretty close to 500, and the V_PulseLoc2 and V_PulseLoc1 values are pretty close to what one gets by hand. This might not be the best way to do this, but it might be helpful.
September 30, 2008 at 08:08 pm - Permalink
October 1, 2008 at 09:01 am - Permalink
Function GradCal(w, Size, Points1D, Points2D, Threshhold, Box, grad)<br />
wave w; variable Size, Points1D, Points2D, Threshhold, Box; string grad; <br />
variable i<br />
string GradLeft=grad+"Left", GradRight=grad+"Right", GradCalib=grad+"Calib"<br />
string GradBW=grad+"BW", GradCent=grad+"Cent"<br />
Make /O/N=(Points2D)/D $GradLeft, $GradRight, $GradCalib, $GradBW, $GradCent<br />
Wave Left=$(GradLeft), Right=$(GradRight), Calib=$(GradCalib), BW=$(GradBW), Cent=$(GradCent);<br />
<br />
for(i=0;i<Points2D;i+=1)<br />
wavestats /R=[Points1D*i,Points1D*(i+1)-1] /Q w<br />
pulsestats/b=(Box) /f=(Threshhold) /L=(V_max, mean(w,Points1D*i,Points1D/2+Points1D*i)) /R=[Points1D*i,Points1D*(i+1)-1] /M=20 /Q w<br />
// printf "%g, %g\r" V_PulseLoc1, V_PulseLoc2<br />
if(i<Points2D/2-1)<br />
Left[i]=V_PulseLoc1; Right[i]=V_PulseLoc2<br />
else<br />
Left[i]=V_PulseLoc2; Right[i]=V_PulseLoc1<br />
endif<br />
endfor<br />
<br />
BW=(Left-Right)*31250/Points1D<br />
Calib=BW/(42.5775*Size)<br />
<br />
End<br />
October 1, 2008 at 02:11 pm - Permalink
Stephen R. Chinn
October 2, 2008 at 03:39 am - Permalink
October 2, 2008 at 07:05 am - Permalink
Yes, Edgestats is pretty much PulseStats for one edge.
Regarding FindPeak, it's rather primitive and needs an upgrade. From the help:
What it lacks are estimates of the expected width, like PulseStats' /M=dx, and amplitude, like FindLevels.
Software Engineer, WaveMetrics, Inc.
October 2, 2008 at 11:27 am - Permalink