Thompson - Computing for Scientists and Engineers (523188), страница 72
Текст из файла (страница 72)
nYou may wish to investigate in more detail the effect of binning on a Lorentzianspectrum.Exercise 10.19Consider the effect of binning a Lorentzian-shaped spectrum extending frominto N bins, which are sometimes called channels in the parlance of multichannel analyzers.(a) Show that the fractional error from binning, E (x), defined as (convoluted —Lorentzian) / (Lorentzian) is given by(10.41)(b) Use the Taylor expansion of the arctangent function to approximate this as(10.42)10.3CONVOLUTIONS AND FOURIER TRANSFORMS397(c) Assume that the binning is done into N = 128 bins, about the lowest resolution that one would typically use. Show that the fractional error at the peak ofthe Lorentzian is only -1.3 × 10-3. nThe increasing width with decreasing height of all the convolutions that we havemade explicitly suggests that convolution preserves area. This is indeed so, as canbe shown from the definition by integrating both sides with respect to x, interchanging the order of integrations on the right-hand side, then changing the variable of integration for y2, keeping x' fixed, to prove the area-conserving property of a convolution, namelyExercise 10.20Show in detail the steps indicated for proving this area-conserving property ofconvolutions.
nAnother observation that you may have made is that, starting with two functionsthat are symmetric about the origin, such as the Lorentzian (Figure 10.5) and theboxcar (Figure 10.8), their convolution is also symmetric about the origin. Is this ageneral property for convolutions? Try it and see.Exercise 10.21Suppose that y1 and y2 are symmetric about the origin, yi (-x) = yi (x) fori = 1,2. Use this property in the convolution definition (10.34), together with achange of variable of integration from x' to -x', to show that their convolutionis also symmetric about the origin. nThe filtering and windowing of Fourier integral transforms, of which the boxcarconvolution provides an example, proceeds similarly to that discussed in Section 9.8.
The techniques are especially important in the reconstruction of three-dimensional images by tomography, as discussed at an introductory level by Solomonin his teaching module. The subject of filtering is covered in, for example, Hamming’s book on digital filters. Bracewell’s text on Fourier transforms and their applications makes extensive use of convolutions.We now have some insight into the properties of convolutions and into possibleapplications of them. There are few interesting examples of analytical functionswhose convolutions are straightforward to calculate from the definition (10.34).Therefore, we proceed to two alternatives; first, numerical estimation of convolutions of discretized functions, then convolutions performed analytically by usingFourier integral transforms.398FOURIER INTEGRAL TRANSFORMSProgram for convoluting discretized functionsIn most practical applications the functions that one wants to convolute are experimental data, necessarily at discrete points.
Therefore, the integral convolutions thatwe have discussed so far are idealizations when applied to data. It is therefore useful, and even interesting, to construct the analogous discretized convolution. Thiswill bear a similar relation to the integral definition of the convolution, (10.34), asthe discrete Fourier transform in Section 9.2 bears to the Fourier integral transformin this chapter.We write the discretized convolution as(10.44)where array notation is used both for the two inputs and for the output convolution.There is an ambiguity in such a definition, because array elements less than zero (inC, or less than unity in Fortran) are often problematic.
In the program below, Imade the simplest, and often realistic, choice of ignoring terms which have the indexvalue [ i - j ] < 1. This is referred to as “truncation.” An alternative choice,called “wraparound,” is to associate with zero or negative index values the valuesread down from the highest index in the array of y2 values.The discrete convolution function convolve_arrays given below implements(10.44) assuming truncation. In order to test the function, the main program waswritten to allow preparation of input arrays, and also output to a file for graphing theresults. The file CONARRAY contains the N values of y1, and CONARRAY2 contains the values of y2. The discrete convolution is output to the console and to thefile CONV.
The range of i values for which the convolution in (10.44) is computedis from imin to imax, inclusive. No checking is done to ensure that this range issensible.PROGRAM 10.1 Direct discretized convolution of two arrays according to (10.44).#include <stdio.h>#include <math.h>#define MAX 513main()/* Convolute Arrays *//* Input arrays then convolute them */FILE *fin,*fout,*fin1,*fin2;FILE *fopen();doubley1[MAX],y2[MAX],y12[MAX],yin,yin1,yin2;int N,i,j,imin,imax;10.3 CONVOLUTIONS AND FOURIER TRANSFORMSchar yn,ynfile,wa;void convolve_arrays();printf("Convolute Arrays: Input # of data, N:\n");scanf("%i",&N);if( N > MAX-1 )printf ("!! N=%i exceeds array sizes %i\n",N,MAX-1);exit(l) ;}printf("Prepare input file(s)? (y or n): ");scanf("%s",&yn) ;if ( yn == 'y' )for ( i = 1; i <= 2; i++ ){printf("\nPrepare y%i file? (y or n): ",i);scanf("%s",&ynfile);if ( ynfile == 'y' )switch (i) /* Open file for write */case 1: fout = fopen("CONARRAY1","w"); break;case 2: fout = fopen("CONARRAY2","w"); break;printf("Input %i data for file %i:\n",N,i);for ( j = 1; j <= N; j++ )printf ("%i: ",j):scanf("%lf",&yin) ;switch (i)case 1: y1[j] = yin;case 2: y2[j] = yin;}fprintf(fout,"%lf\n",yin);/* Ready to reuse */}fclose(fout); rewind(fout) ;}}} /* end prepare input files */printf("\nConvolutioncalculation\n");printf("Write over output (w) or Add on (a): ");scanf("%s",&wa);fout = fopen("CONV",&wa);399400FOURIER INTEGRAL TRANSFORMSfin1 =fopen("CONARRAY1", "r") ;fin2 = fopen("CONARRAY2", "r");for ( j = 1; j <= N; j++ )/* Data input from arrays */{fscanf(fin1,"%lf\n",&yin1); y1[j] = yinl;fscanf(fin2, "%lf\n",&yin2); y2[j] = yin2;}printf("Choose min (>=l) & max (>=min,<=N) index:\n");scanf ("%i%i" ,&imin,&imax);/* Convolution */convolve_arrays(N,imin,imax,yl,y2,y12);for ( i = imin; i <= imax; i++ )/* Output convolution *i{printf("\n%i %8.4lf",i,y12[i]);fprintf(fout,"%i %lf\n",i,yl2[i]);}printf("\n\nConvolution of arrays ends");voidconvolve_arrays(N,imin,imax,yl,y2,y12)/* Convolute arrays yl & y2 for index imin to imax */Array indices <= 0 do not contribute *//*double y1[],y2[],y12[];int N,imin,imax;{int i,j;for ( i = imin; i <= imax; i++ )yl2[i] = 0;for ( j = 1; j <= N; j++ ) /* Accumulate convolution */if ( i > j ) yl2[i] = yl2[ij+y1[j]*y2[i-j];This convolution program, Convolute Arrays, appears to be mostly a mainprogram for file handling.
The program will, however, be useful for performingdiscretized convolutions of Gaussians and Lorentzians (whose integral transformswe consider below) and for convoluting data that you may have.Exercise 10.22In order to test convolve_arrays and the main program, Convolute Arrays,one may check that the dependence of an integral convolution, (10.34). on x10.3 CONVOLUTIONS AND FOURIER TRANSFORMS401and of a discretized convolution, (10.44), on i are similar when the steps in xare small compared to the rate of variation of the functions y1 and y2.To do this, for y2 use the Lorentzian function withand x in the rangefrom -3.5 to +3.5 by steps ofthis has N = 36.
For y1 use the boxcar function, Figure 10.8, with L = 1. After making use of the correspondencex = -3.5 + (i - 1)verify that the shape of this curve (but not the verticalscale) is the same as the curve for L = 1 in Figure 10.9 for the integral convolution, except near the endpoints where the truncation assumption spoils theagreement. nAn advanced treatment of convolution algorithms that is suitable for use with theFFT is given in Nussbaumer’s book.Now that we have some experience with numerical convolutions and insight intotheir properties, let us return to the analytical developments.Fourier integral transforms and convolutionsWe now derive analytical results for convolutions that depend on properties of Fourier integral transforms.
These results will provide a powerful means of calculatingconvolutions, as we will show for Gaussian and Lorentzian functions.The general result is the Fourier corvolution theorem. This theorem states thatthe Fourier integral transform of the convolution of two functions is proportional tothe product of their Fourier transforms. In symbols the convolution theorem reads(10.45)where the integral(10.46)denotes the Fourier integral transform of the convolution of the two functions,whose Fourier transforms are Y1(k) and Y2(k),Exercise 10.23To prove the convolution theorem, first write down Y1(k) and Y2(k) in termsof their definition, (10.4), using integration variables x' and x". Then changethe x" variable to x - x' and perform the x' integration first, which producesthe convolution of y1 with y2.
Show that the integral that remains is just proportional to the Fourier transform of this convolution, as (10.46) claims. nThis remarkable theorem has an immediate use for calculating convolutions. Byexpressing the convolution in terms of its Fourier transform using (10.45), then using the convolution theorem (10.46) for the transform you will find immediately that402FOURIER INTEGRAL TRANSFORMS(10.47)Exercise 10.24Find this result for a convolution in terms of the Fourier integral transforms ofthe two functions by carrying through the indicated steps. nWe immediately put this convolution formula to use in computing the convolution of two Gaussians or of two Lorentzians.Convolutions of Gaussians and of LorentziansYou probably wondered why, when working examples of convolutions at the beginning of this section, we did not use the common Gaussian and Lorentzian functions,of which I seem to be so fond.
Although convolution of Gaussians can be done directly from the definition (10.34), the derivation is messy. Convolution of two Lorentzians can be carried out directly, but it is also tedious unless complex-variablemethods are used.Convolution of two Gaussians by use of the convolution theorem, (10.46), isvery direct because the Fourier transform of a Gaussian is just a Gaussian, as weshowed in Section 10.2.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
