PEAKFIND (Feb93) grasp.grtools PEAKFIND (Feb93) NAME peakfind -- Fit peaks in helioseismic power spectra and L-Nu diagrams. USAGE peakfind input output table PARAMETERS input Input power spectrum image. output Output fit image. table Output Mode-frequency table. lm_start = INDEF Starting L or M value to be fit. lm_end = INDEF Ending L or M value to be fit. n_start = 1 Starting N-value to be fit. n_end = 10 Ending N-value to be fit. window = "" Window function power spectrum. fit_background = yes Fit the background. background = "" Background image name if not fitting. restart = 1 Number of simplex restarts per fit. nleaks = 0 Number of leaks to include in the fit. DESCRIPTION PEAKFIND uses a nonlinear, maximum likelihood technique to fit lorentz profiles to peaks in helioseismic power spectra (M-Nu or L-Nu diagrams). For each ridge in the range lm_start => lm_end and n_start => n_end, the task determines the mode frequency, FWHM and amplitude, of primary mode plus the closest nleaks leakage modes. If fit_background = yes, a linear background is also fit to the ridge. If fit_background = no, a predetermined background image, bkg_spectrum may be used. The input power spectrum is expected to be in GONG format. The task GRTOOLS.SETMWCS may be used to put the required MWCS header into an image that was not generated by GRASP. The type of power spectrum is given by the input header keyword DTYPE and may be either LNU or PS for L-Nu diagrams and M-Nu diagrams respectively. The output image is a 3-d image with band one containing the derived limit spectrum, and band 2 containing the original input power spectrum. Table is a STSDAS table containing the results of the fit. The table may be operated upon by any of the tasks of the TABLES package. The table columns are: Column Data Type Default Fmt Plot Label Comment n I %2d "" n-value l I %4d "" l-value m I %4d "" m-value leak I %2d "" 0=primary mode; 1=leak nu R %10.4f microHz frequency dnu R %10.4f microHz frequency error fwhm R %8.4f microHz FWHM dpfwhm R %8.4f microHz FWHM positive error bar dmfwhm R %8.4f microHz FWHM negative error bar psamp R %10.4f (cm/s)**2 Power spectrum amplitude dppsamp R %10.4f (cm/s)**2 " " " positive error bar dmpsamp R %10.4f (cm/s)**2 " " " negative error bar tsamp R %10.4f cm/s Mode amp in time domain dptsamp R %10.4f cm/s " " " positive error bar dmtsamp R %10.4f cm/s " " " negative error bar bkg0 R %10.4f (cm/s)**2 Background intercept dbkg0 R %10.4f (cm/s)**2 " " error bkg1 R %10.4f (cm/s)**2 Background slope dbkg1 R %10.4f (cm/s)**2 " " error merit D %8 "" Merit Function The errors on the FWHM and Amplitude/Power are evaluated with respect to the logarithm of those parameters. The quoted "+/-" error is with respect to the log distribution. The error is then translated to the linear distribution resulting in separate (i.e., unequal) "+" and "-" error bars. Transformation of power spectrum amplitudes to time domain is not yet enabled. If window = "", the a perfect observing window is assumed. However, window may be given the name of the observing window power spectrum image created by the task PWINDOW, in which case the non-perfect window is taken into account while fitting the modes. The maximum likelihood function is minimized using the Downhill Simplex method. The number of simplex restarts is given by restart. The fitting range extends to half the distance (in frequency) to the next ridge. The initial guesses for the mode parameters are determined as follows: Guess frequencies: The guess frequencies are determined using spline interpolation in a sparse table of the South Pole 1987 results. Guess FWHM The guess FWHM are computed from power-law graphical fits to Libbrecht's curve as a function of frequency. Dependence on l taken from linear fit to Jefferies et al. The answer constrained to be at least 0.1 Guess Amplitudes The amplitudes are determined from the mean (M-Nu) or maximum (L-Nu) value of the data within 3 pixels of the guess frequency. Guess Background intercept and slope The guess background parameters are determined from a linear fit to the mean value of the two groups of 5 pixels at the ends of the fitting range. EXAMPLES 1. To fit all the peaks and background for 5 <= n <=10 and 50 <= l <= 100 using the window function "pwindow": gr> peakfind inimage outimage outtable lm_s=50 lm_e=100 n_s=5 n_e=10 TIME REQUIREMENTS On a SPARCstation 10, the average time to fit a ridge including 8 leaks and the background with one simplex restart is 9.5 cpu seconds for the South Pole 1987 data set (1 month resolution plus a window function). BUGS No phase information is currently used from the input or saved to the output. This version of PEAKFIND only takes into account L-leakage. Thus, it cannot yet be used to fit low-L modes (l <= 20) which suffer from n-leakage as well. This version of PEAKFIND is not smart enough to determine that a peak is "not there". Thus, it will result in a fit to noise spikes where peaks do not really exist. This version of PEAKFIND has problems with unresolved leaks (i.e., fitting high frequency and high L modes). This version of PEAKFIND has not been tested on M-Nu diagrams due to a lack of GONG data of sufficient length in time. Any suggestions from the community on how to deal with any of the above problems are most welcome. SEE ALSO pwindow, setmwcs "Modeling, of, Solar, Oscillation, Power, Spectra" Anderson, E.R., Duvall, T.L., Jr., and, Jefferies, S.M. Ap.J., December, 1, 1990, (IN, PRESS).