; docformat = 'rst'
;+
; :Author: Paulo Penteado (pp.penteado@gmail.com), Feb/2010
;-
;+
; :Description:
; Examples to test and illustrate the different behaviors of pp_convol.
;
; Makes up a high resolution spectrum and convolves it to lower resolutions.
;
; :Params:
; width : in, optional, default=0.04
; The width for the convolution kernels.
;
; :Examples:
; Convolve spectra to different resolutions::
;
; .com pp_convol ;pp_convol_test is in pp_convol.pro, so it will not be found by itself
; pp_convol_test,0.04
; pp_convol_test,0.02
;
; .. image:: pp_convol_test_1.png
; .. image:: pp_convol_test_2.png
;
; :Uses: pp_resample,pp_integral
;
; :Author: Paulo Penteado (pp.penteado@gmail.com), Feb/2010
;-
pro pp_convol_test,width
;Example and test case for pp_convol
compile_opt idl2
;Defaults
width=n_elements(width) eq 1 ? width : 0.04
;Make up a high resolution spectrum and plot it
nx=1001
x=dindgen(nx)/(nx-1d0)
noise=0.005
y=1d0+randomu(seed,nx,/normal)*noise
nlines=30
df=0.7d0
depths=(randomu(seed,nlines)-0.5d0)*df
wf=0.025
widths=randomu(seed,nlines)*wf
centers=randomu(seed,nlines)
for i=0,nlines-1 do y+=depths[i]*exp(-((x-centers[i])/widths[i])^2)
iplot,x,y,name='Original spectrum',/insert_legend,$
title='Original at '+strtrim(string(x[1]-x[0]),2)+', convolved to '+strtrim(string(width),2)
;Convolve with Gaussians at lower resolution
yc0=pp_convol(y,x,/gaussian,width=width,/local)
iplot,x,yc0,/over,color=[255,0,0],name='Gaussian (literal)',/insert_legend
yc0=pp_convol(y,x,/gaussian,width=width)
iplot,x,yc0,/over,color=[0,0,255],name='Gaussian (sampled)',/insert_legend
yc0=pp_convol(y,x,/step,width=width,/local)
;Convolve with rectangular kernels at lower resolution
iplot,x,yc0,/over,color=[0,255,0],name='Step (literal)',/insert_legend
yc0=pp_convol(y,x,/step,width=width)
iplot,x,yc0,/over,color=[255,0,255],name='Step (sampled)',/insert_legend
end
;+
; :Description:
; Convolves the provided y(x) function with either a rectangular or Gaussian
; kernel of the given width, or the provided arbitrary kernel (not yet implemented).
;
; The x grid does not need to be regular, and, contrary to common practice, no resampling
; is done to a regular grid to calculate the convolution (resampling is not needed). Also
; contrary to common practice, for rectangular or Gaussian kernels the integral is done
; analytically, not numerically. In the case of a Gaussian kernel, it extends from the Gaussian
; center up to the point where the terms in the integral become 0 (at double precision).
;
; Also, this routine allows for both literal and sampled interpretations of the input function (set
; with the local keyword).
;
; :Returns:
; The result of the convultion is given for the same x locations as the input,
; there is no automatic resampling. If resampling is intended after the convolution,
; pp_resampled is recommended, as it preserves the areas.
;
; :Params:
; x : in, required
; An array of locations where the function is sampled. Must be ordered (increasing or decreasing).
; y : in, required
; An array with the function values corresponding to the locations in x.
;
; :Keywords:
; step : in, optional, default=1
; If set, the convolution kernel is a rectangular function of the given width.
; gaussian : in, optional, default=0
; If set, the convolution kernel is a Gaussian, of fwhm of the given width.
; width : in, required
; The width of the convolution kernel (for step or gaussian).
; grid_tolerance : in, optional,default=1d-6
; Not yet implemented. Specifies the maximum fractional variation in the input grid to accept it as uniform.
;
; If the input grid is found to be uniform, the integration can be done faster.
; kernelx : in, optional
; Not yet implemented. An array of the x locations (centered at x=0) of the arbitrary convolution kernel to use.
;
; Must be ordered in increasing x.
; kernely : in, optional
; Not yet implemented. An array of the the arbitrary convolution kernel's values to use.
; Must make a function of unit area.
; local : in, optional, default=0
; If not set, the y values are interpreted as a measured sample, that is, as the average of the "flux" falling inside the region
; centered at the corresponding x values: x locations are interpreted as the centers of bins. Since this interpretation is equivalent
; to a literal interpretation of a function made of rectangular regions, in this method there are no interpolations (partial areas, if
; any, are just the corresponding fractions of the rectangles).
;
; If set, interprets the input function literally.
; newton : in, optional, default=0
; Passed on to pp_integral. If local is set and newton is set, the function integrations are done with int_tabulated, which
; uses a 5 point Newton-Cotes formula. If local is not set, this has no effect.
; lsquadratic : in, optional, default=0
; Passed on to pp_integral. If the method requires interpolation, this is passed on to interpol,
; to select 4 point quadratic interpolation.
;
; If none of lsquadratic, quadratic and spline are set, interpolation is linear.
; quadratic : in, optional, default=0
; Passed on to pp_integral. If the method requires interpolation, this is passed on to interpol,
; to select 3 point quadratic interpolation.
;
; If none of lsquadratic, quadratic and spline are set, interpolation is linear.
; spline : in, optional, default=0
; Passed on to pp_integral. If the method requires interpolation, this is passed on to interpol,
; to select spline interpolation.
;
; If none of lsquadratic, quadratic and spline are set, interpolation is linear.
; xinc : out, optional
; Returns the x values, in increasing order.
; yinc : out, optional
; Returns the y values, in order of increasing x.
; xstart : out, optional
; If local is not set, returns the location where each input pixel starts. If local is set, this is ignored.
; xend : out, optional
; If local is not set, returns the location where each input pixel ends. If local is set, this is ignored.
; xbox : out, optional
; Used when not in local mode, to return the input x values in pairs of equal values, so that
; plotting xbox,ybox shows the rectangles that are how the input values were interpreted.
; ybox : out, optional
; Used when not in local mode, to return the input y values in pairs of equal values, so that
; plotting xbox,ybox shows the rectangles that are how the input values were interpreted.
;
; :Examples:
;
; Make a well-sampled input function y(x) with fwhm=1::
;
; x=(12d0*dindgen(1001)/1d3)-6d0
; y=exp(-(x^2*4d0*alog(2d0)))*8d0
; iplot,x,y,/isotropic,name='Original function of fwhm=1',/insert_legend,thick=2.
;
; Do the convolution and look at the results::
;
; yc1=pp_convol(y,x,/step,width=4d0)
; iplot,x,yc1,/over,color=[255,0,0],name='Convolution width=4, rectangular kernel (sampled)',/insert_legend
; yc2=pp_convol(y,x,/step,/local,width=4d0)
; iplot,x,yc2,/over,color=[0,0,255],name='Convolution width=4, rectangular kernel (literal)',/insert_legend
; yc3=pp_convol(y,x,/gaussian,width=4d0)
; iplot,x,yc3,/over,color=[0,255,0],name='Convolution width=4, Gaussian kernel (sampled)',/insert_legend
; yc4=pp_convol(y,x,/gaussian,width=4d0,/local)
; iplot,x,yc4,/over,color=[255,0,255],name='Convolution width=4, Gaussian kernel (literal)',/insert_legend
;
; .. image:: pp_convol.png
;
; See also the example in pp_convol_test, for a comparison of the methods with a more realistic (spectrum-like) function.
;
; :Todo: Implement arbitrary kernel, grid_tolerance, and quadratic and cubic Gaussian integrations.
;
; :Uses: pp_resample,pp_integral
;
; :Author: Paulo Penteado (pp.penteado@gmail.com), Feb/2010
;-
function pp_convol,y,x,step=step,gaussian=gaus,width=width,grid_tolerance=gtol,$
kernelx=kx,kernely=ky,$
local=local,newton=newton,$
lsquadratic=lsquadratic,quadratic=quadratic,spline=spline,$
xinc=ix,yinc=iy,xstart=xstart,xend=xend,xbox=xbox,ybox=ybox
compile_opt idl2
;Defaults
local=n_elements(local) eq 1 ? local : 0
step=n_elements(step) eq 1 ? step : 1
gaus=n_elements(gaus) eq 1 ? gaus : 0
nkx=n_elements(kx) & nky=n_elements(ky)
usekernel=((nkx eq nky) && (nkx gt 0L)) ? 1 : 0
if (gaus) then begin
step=0
usekernel=0
endif
if (usekernel) then begin
step=0
gaus=0
endif
gtol=n_elements(gtol) eq 1 ? gtol : 1d-6
;Check that the dimensions of x and y match and are acceptable
nx=n_elements(x) & ny=n_elements(y)
if (nx ne ny) || (nx eq 0) then message,'x and y provided must be arrays of the same length'
;Make ix and iy to work with, 1D and in increasing order (x is assumed sorted, but can be increasing or decreasing order)
if (x[nx-1] lt x[0]) then begin
ix=reverse(reform(x)) & iy=reverse(reform(y))
endif else begin
ix=reform(x) & iy=reform(y)
endelse
;Determine whether the input grid is regular
dmax=max(ix[1:*]-ix[0:nx-2],/abs,min=dmin)
regular=((dmax-dmin)/dmin le gtol)
if (step) then begin ;Use a constant function as the kernel
;Calculate the locations of the start and end of the interval that will go into each output point
xmin=x-width/2d0
xmax=x+width/2d0
;Calculate the input function's integral over the proper regions
areas=pp_integral(x,y,xmin=xmin,xmax=xmax,newton=newton,local=local,$
lsquadratic=lsquadratic,quadratic=quadratic,spline=spline)
ret=areas/width
endif
if (gaus) then begin ;Use a Gaussian function as the kernel
maxwid1=27.3d0 ;In double precision, exp(-x^2) underflows between x=27.29d0 and x=27.3d0
maxwid0=8.27d0 ;In double precision, gauss_pdf(x) saturates (becomes 1d0) between x=8.26d0 and x=8.27d0
;maxwid0=5.888d0 ;In double precision, erf(x) saturates (becomes 1d0) between x=5.8870 and x=5.888d0
a=width^2/(4d0*alog(2d0)) ;The constant that will go into the exponential (exp(-x^2/a))
g=1d0/sqrt(a*!dpi) ;Constant that multiplies the exponential in the Gaussian (g*exp(-x^2/a))
wid0=sqrt(a/2d0)*maxwid0
wid1=sqrt(a/2d0)*maxwid1
;Calculate the integral
;This is not numerical integration: it is analytical integration of a function made of a finite number of pieces
if (local) then begin ;Literal interpretation of the function
;Calculate the locations of the start and end of the interval that will go into each output point
minloc1=(value_locate(ix,ix-wid1))>0L
maxloc1=(value_locate(ix,ix+wid1)+1L)<(nx-2L)
minloc0=(value_locate(ix,ix-wid0))>0L
maxloc0=(value_locate(ix,ix+wid0)+1L)<(nx-2L)
;The integral is of y(x)*g(x-x0)*dx, with y(x)=c0+c1*x
;Terms to calculate c0 and c1, coefficients of y(x)=c0+c1*x
y1=iy[1:nx-1] & y0=iy[0:nx-2] & x1=ix[1:nx-1] & x0=ix[0:nx-2]
c0=(y0*x1-y1*x0)/(x1-x0)
c1=(y1-y0)/(x1-x0)
;Integral of the first term (c0*g(x-x0)*dx)
w=where((maxloc0 gt minloc0),nw)
ret0=dblarr(nx)
for i=0,nw-1 do begin
j=w[i]
xl=ix[j]
rng=lindgen(maxloc0[j]-minloc0[j]+1)+minloc0[j]
x0l=(x0[rng]-xl)/sqrt(a)
x1l=(x1[rng]-xl)/sqrt(a)
d=(c0+c1*xl)[rng]
ret0[j]=total(d*(erf(x1l)-erf(x0l)))
;ret0[j]+=total(d*(gauss_pdf(x1l)-gauss_pdf(x0l)))
endfor
;Integral of the second term (c1*x*g(x-x0)*dx)
w=where((maxloc1 gt minloc1),nw)
ret1=dblarr(nx)
for i=0,nw-1 do begin
j=w[i]
xl=ix[j]
rng=lindgen(maxloc1[j]-minloc1[j]+1)+minloc1[j]
x0l=(x0[rng]-xl)/sqrt(a)
x1l=(x1[rng]-xl)/sqrt(a)
ret1[j]+=total(c1[rng]*(exp(-x0l^2)-exp(-x1l^2)))
endfor
ret=g*0.5d0*(sqrt(!dpi*a)*ret0+a*ret1)
;ret=g*(ret0+0.5d0*a*ret1)
endif else begin
;Calculate the locations of the start and end of the interval that will go into each output point
xstart=[ix[0]-(ix[1]-ix[0])/2d0,(ix[0:nx-2]+ix[1:nx-1])/2d0]
xend=[xstart[1:nx-1],ix[nx-1]+(ix[nx-1]-ix[nx-2])/2d0]
minloc0=(value_locate(xstart,xstart-wid0))>0L
maxloc0=(value_locate(xend,xend+wid0)+1L)<(nx-1L)
;The integral is of y(x)*g(x-x0)*dx, with y(x)=c0
xsig=ix/sqrt(a/2d0) ;x normalized by the Gaussian standard deviation
xstartsig=xstart/sqrt(a/2d0) ;x normalized by the Gaussian standard deviation
xendsig=xend/sqrt(a/2d0) ;x normalized by the Gaussian standard deviation
w=where((maxloc0 gt minloc0),nw)
ret=dblarr(nx)
for i=0L,nw-1 do begin
j=w[i]
ret[j]+=total(iy[minloc0[j]:maxloc0[j]]*(gauss_pdf(xendsig[minloc0[j]:maxloc0[j]]-xsig[j])-gauss_pdf(xstartsig[minloc0[j]:maxloc0[j]]-xsig[j])))
endfor
endelse
endif
if (~local) then begin
;Calculate the start and en of each output point
xstart=[ix[0]-(ix[1]-ix[0])/2d0,(ix[0:nx-2]+ix[1:nx-1])/2d0]
xend=[xstart[1:nx-1],ix[nx-1]+(ix[nx-1]-ix[nx-2])/2d0]
;Create xbox and ybox arrays, with constant y values spanning the start and end of each x point
if (arg_present(xbox) || arg_present(ybox)) then begin
nix=n_elements(ix)
xbox=dblarr(2L*nix) & ybox=dblarr(2L*nix)
xbox[0:*:2]=xstart & xbox[1:*:2]=xend
ybox[0:*:2]=iy & ybox[1:*:2]=iy
endif
endif
return,ret
end