Thompson - Computing for Scientists and Engineers (523188), страница 66
Текст из файла (страница 66)
Program speed could be increased byusing machine language to perform the bit reversal.The Fast Fourier Transform program is self-testing by use of the symmetry between the discrete Fourier transform relation (9.6) and the inverse relation(9.7). All that needs to be done is to reverse the sign of the variable IT, which hasIT = + 1 to obtain the coefficients ck from the data yj, and IT = -1 to invert the ck9.7PROGRAM FOR THE FAST FOURIER TRANSFORM361and obtain the yj Recall that, according to Sections 9.2 and 9.3, our derivation assumes that the data are real variables.PROGRAM 9.3 Radix-2 fast Fourier transform checking program.#include <stdio.h>#include <math.h>#define MAX 33main(){/* Fast Fourier Transform;radix 2; Checking program */double yreal [MAX] ,yimag [MAX] ;/* Fortran compatibility; arrays are used from [l] */double yin,twopi;int nu,N,MAX1,j,it,iftest;void FFT();twopi = 8*atan(1.O);printf("FFT: Input nu: ");scanf ("%i",&nu) ;N = pow(2,nu); MAX1 = MAX-1;printf("nu = %i, so # data N = %i\n",nu,N);if(N>MAX1)printf("\n !! N=%i > MAX1=%i?n",N,MAX1);printf("End FFT"); exit(l);printf("\nInput %i y data\n",N);for ( j = 1; j <= N; j++ )scanf ("%f",&yin) ;yreal[j] = yin; yimag[j] = 0;it = 1; /* FFT from y back into y *//* Warning; Original y is destroyed; copy if needed */FFT(nu,N,it,twopi,yreal,yimag);printf("\nOutput of FFT; yreal, yimag\n");for ( j = 1; j <= N; j++ ){printf("\n%g%g",yreal[j],yimag[j]);}it = -1; /* Inverse FFT from y back into yshould approximately restore y */362DISCRETE FOURIER TRANSFORMS AND FOURIER SERIESFFT(nu,N,it,twopi,yreal,yimag);printf("\nOutput of FFT; yreal, yimag\n");printf("{yimag should be zero to machine tolerance\n}");for ( j = 1; j <= N; j++ )printf("\n%g%g",yreal[j],yimag[j]);iftest = 1; /* FFT speed test option; manual timing */while ( iftest == 1 )printf("\nInput a 1 for FFT speed test: ");scanf("%i",&iftest);if ( iftest == 1 )printf("Input nu, < 11 : "); scanf("%i",&nu);N = pow(2,nu);printf("Start timing now\n");for ( j = 1; j <= N; j++ )yreal[j] = j; yimag[j] = 0;it = 1; /* FFT */FFT(nu,N,it,twopi,yreal,yimag);printf("Stop timing now\n");elseprintf("End FFT"); exit(0);} /* end iftest speed loop */voidFFT(nu,N,it,twopi,yreal,yimag)/* Fast Fourier Transform of N data in yreal+i(yimag) */double yreal[],yimag[];double twopi;int nu,N,it;doubletpn,angle,wr,wi,tempr,tempi,norm;intii,inu,nu1,nsub,i,j,inup1,in1,jsort;int bitrev();tpn = twopi/N;ii = -it;inu=0; nu1=nu-1; nsub=N/2; /* set up first subinterval *//* Loop over levels; lengths N/2, N/4, .
. . */9.7PROGRAM FOR THE FAST FOURIER TRANSFORMfor ( i = 1; i <= nu; i++ ){while ( inu < N )for ( j = 1; j <= nsub; j++ ){ /* Transform subinterval */angle = tpn*bitrev((inu/pow(2,nu1)),nu);wr = cos (angle) ; wi = ii*sin(angle);inup1 = inu+1; in1 = inup1+nsub;tempr = yreal[in1]*wr - yimag[in1]*wi;tempi = yreal[in1]*wi + yimag[in1]*wr;yreal[in1] = yreal[inup1]-tempr;yimag[in1] = yimag[inup1] -tempi;yreal[inup1] = yreal[inup1]+tempr;yimag[inup1] = yimag[inup1]+tempi;inu = inu+l;inu = inu+nsub;inu=0;nu1=nu-1;nsub=nsub/2; /* next subinterval *//* Reverse bit pattern and sort transform */for ( j = 1; j <= N; j++ )jsort = bitrev(j-l,nu)+ 1;if ( jsort > j )tempr =yreal[j]tempi =yimag[j]yreal[j]; /* swap values */= yreal[jsort]; yreal[jsort] = tempr;yimag[j];= yimag[jsort]; yimag[jsort] = tempi;norm = 1/sqrt( (double) N); /* Normalize */for ( j = 1; j <= N; j++ ){yreal[j] = norm*yreal[j];yimag[j] = norm*yimag[j];int bitrev(j,nu)/* Bit reverse the nu bits of j into jr */int j,nu;363364DISCRETE FOURIER TRANSFORMS AND FOURIER SERIESint jt,i,jd,jr;jt = j; jr = 0;for ( i = 1; i <= nu; i++ )jd = jt/2;jt = jd;jr = 2*jr + jt-2*jd;return jr;The following exercise uses the analytically-calculated discrete transform for theharmonic oscillator (single frequency) derived in Section 9.2.
It illustrates the usefulness of analytical special cases for checking numerics and programs.Exercise 9.34(a) RunBrogram9.1, Discrete Fourier Transform for Oscillator,from Section 9.2 for N = 16, 64, and 512, with your choice of the oscillatorangular frequency k0, in order to generate the discrete Fourier transform coefficients ck for k = 0 to N - 1. Save these coefficients in a file.(b) Run Fast Fourier Transform to compute the inverse Fourier transformof the array of ck values generated in (a). That is, input IT = -1 and get thedata from the file.
You should recover a complex exponential having angularfrequency k0. Check the accuracy of this inversion. nAnother way to check the FFT program is to generate the discrete Fourier coefficients for the decaying exponential (Section 9.2), which is quite different from thatin Exercise 9.34, thus providing an alternative check of the FFT program.Exercise 9.35(a) Convert (9.20) into expressions for the real and imaginary parts of the ck forexponential decay.
Then code and spot check the program. For example, youcan check the symmetry condition (9.10). Use N = 16, 64, and 512, with yourchoice of the stepsize h, in order to generate the discrete Fourier transform coefficients ck for k = 0 to N - 1. Save these coefficients in a file.(b) Run the program Fast Fourier Transform to compute the inverseFourier transform of the array of ck values generated in (a). That is, inputIT = -1 and get the data from the file created by the program for exponentialdecay. You should recover a decaying exponential at the stepsize h. By takingthe natural logarithm of the resulting coefficients you can check how well the inverse transform perfoims. nSpeed testing the FFT algorithmSpeed testing of the FFT algorithm in order to verify the scaling rule indicated in(9.37) is given as a worked example in Exercise 9.12.
The program Fast F o u r9.8FOURIER ANALYSIS OF AN ELECTROENCEPHALOGRAM365ier Trans form has timing as an option in the main program after the self test.You should now run the timing tests yourself.Exercise 9.36Use the timing (done manually) of the FFT algorithm coded in the programFast Fourier Transform. Use N = 128,256, and 1024 if you are using apersonal computer, Values about factors of 4 or 8 larger than these may be appropriate for a workstation or for a desktop computer with an accelerator board.(If necessary, increase the maximum array size, MAX, to accommodate the largervalues of N.) From your timing results estimate a range for Ta in (9.37) andthus estimate the FFT time for a data sample of N = 8192 = 213 points. nWith the experience gained from running Fast Fourier Transform youshould be ready to analyze real data that do not have an obvious pattern of frequencies.
Such are the data from electroencephalograms.9.8 PROJECT 9B: FOURIER ANALYSIS OFAN ELECTROENCEPHALOGRAMThis project provides an introduction to practical Fourier analysis and it illustrateshow the discrete Fourier transform theory in Section 9.2 is applied to the analysis ofreal data. The electroencephalogram (EEG, or “brainwave”) data are analyzed by theFFT algorithm from Section 9.3 with the program provided in Section 9.7.To the untrained observer the output voltage of the EEG as a function of time indicates very little pattern. The analysis of the EEG data transforms it from the timedomain to the frequency domain, in which the dominance of a few frequencies becomes evident, as you will soon discover. Recall that, as shown in Sections 6.2and 9.1, these transformed amplitudes also provide the best fit to the data (in theleast-squares sense) that can be obtained by using a linear superposition of cosinesand sines.Data from many other physiological rhythms may also be analyzed by theFourier expansion methods in this chapter and Chapter 10.
Among the several examples discussed in Cameron and Skofronick’s book on medical physics are magnetoencephalograms and magnetocardiograms, which use ultrasensitive magnetometersto measure the tiny magnetic fields from the head and heart.Overview of EEGs and the clinical recordSince the 1920s the variation of potential differences between points on the scalp asa function of time have been associated with electrical activity of the brain. The earlymeasurements were subject to much controversy, partly because of the small voltages involved, in the microvolt (µV) range, and particularly because direct interpretation was so subjective.
The advent of computers in the 1950s enabled objectiveanalysis of the frequency components of brain waves, even if not objective interpretation. A typical EEG of an adult human is shown in Figure 9.19.366DISCRETE FOURIER TRANSFORMS AND FOURIER SERIESFIGURE 9.19 Electroencephalogram (EEG) voltages over a l-second interval. The data are connected by solid lines, and the approximation of the data by the first sixteen harmonics is shown dotted. The effects of Lanczos filtering are not visible on this scale.The main frequencies in an adult human EEG are 0.3 to 3.5 Hz (calledwaves, predominant in sleeping adults and awake young children), 8 to 13 Hz (awaves, predominant in awake, relaxed adults), and 18 to 30 Hz waves, appearing sporadically during sleep).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
