Thompson - Computing for Scientists and Engineers (523188), страница 67
Текст из файла (страница 67)
The interpretation of such patterns, especially for diagnosis of neural dysfunction, is discussed at an introductory level by Cameron andSkofronick. Instrumentation for EEG data acquisition is described in the book byCromwell, Weibell, and Pfeiffer. The monograph by Spehlmann provides many examples of EEG traces and their clinical interpretation.For our analysis of EEG data, three EEG pattems of patient MAC (Mental Aptitude Confidential) are provided in Table 9.1. For illustration purposes, we use onlya very short l-second interval, conveniently taken over 32 data points each separatedby 1/31 of a second. Therefore, we may use the radix-2 FFT program developed inSection 9.7, choosing the index v = 5, since 25 = 32. In order to achieve a periodic function, the first and last points of the data have been forced to be the same.The dominant frequencies are characteristic of those of the EEGs of adult humans,as discussed above.
Indeed, you can see a dominant rhythm in Figure 9.19 bynoting a fairly regular crossing of the axis that occurs 16 times in 1 sec, so there hasto be a strong amplitude at about 8 Hz, in the -rhythm section.Three data sets are given in Table 9.1 and displayed in Figure 9.19. As youcan see, they are nearly the same, because they differ only in the amount of noisepresent in the data. I inserted the noise voltages artificially, in such a way that V1has the least noise and V3 has the most noise.9.8FOURIER ANALYSIS OF AN ELECTROENCEPHALOGRAMTABLE 9.1 Data sets for the EEG analysis.ITt(s)V1 (t)V2 (t)V3 (t)(µV)(µV)(µV)13.0010.000013.9312.0720.032220.7719.5619.6430.0645-21.62-12.14-17.1840.0967-19.28-11.45-14.6650.129028.1529.7529.4260.161223.0422.1821.6170.1935-32.87-32.32-33.0780.2258-29.30-38.27-32.4990.258017.6118.8218.51100.290318.6326.6421.16110.3225-14.58-17.42-16.00120.35480.38-8.74-2.7 1130.387 145.2135.9740.30140.419320.9127.5622.94150.4516-30.27-27.19-28.26160.4838-25.33-26.7 1-25.02170.51615.179.146.68180.5483-3.66-7.36-6.22190.5806-28.17-26.13-26.86200.61290.976.474.25210.645 135.6930.3133.00220.67743.561.792.15230.7096-28.18-30.39-29.58240.7419-11.30-8.54-9.2 1250.774 114.2510.6312.94260.80644.56-1.780.39270.8387-0.48-6.38-3.90280.870914.9916.6717.13290.903219.6822.8421.56300.9354-2.171.69-1.71310.9677-16.35-20.65-18.00321.000013.9312.0713.00367368DISCRETE FOURIER TRANSFORMS AND FOURIER SERIESA characteristic of noise in data is that, because it is random from one data pointto the next, it appears to have a frequency about equal to the reciprocal of the sampling interval; thus, here the noise frequency should be about 30 Hz.
You may testthis by seeing whether the voltages reconstructed from the Fourier amplitudes, butomitting the high-frequency components, are in essential agreement with each other.In these EEG data, the average value (the DC level of the EEG voltage) has been removed. Thus, the a0 coefficient in the frequency spectrum should be zero. Nowthat we have the EEG data, we are ready to make Fourier analyses.Program for the EEG analysisIn this subsection we construct a program for analyzing the EEG data in Table 9.1by modifying the FFT program developed in Section 9.7.
Because there are severaldata sets and filtering options to choose from, I have included file input and analysisselections within the program.The program EEG/ FFT consists of three main sections:Optional preparation of EEG data input to file EEGVn, where n = 1, 2, or 3.(i)Analysis of one of these EEG data sets, including the one just just prepared,(ii)to obtain the Fourier amplitudes, ck.
These amplitudes are written to file EEGoutNwhere n = 1, 2, or 3.(iii) Remaking of the EEG with three options for filtering:1. No filtering.2. Lanczos filtering, as explained in the following subsection on filteringthe EEG.3. Truncated filtering, with equal weight up to maxf, then zero weight.For each filtering option, the remade EEG is written to file FilterEEGn, wheren = 1, 2, or 3. The file may either be rewritten or added to. These options allowyou to explore the EEG analysis and to save the results for input to a graphics application, as shown in Figure 9.19. The source program is as follows, with the FFTand bit rev functions being just those given in Section 9.7, so they are omittedfrom the program listing.PROGRAM 9.4 EEG FFT analysis, with variable data and filtering options.#include <stdio.h>#include <math.h>#define MAX 33main()/* EEG FFT analysis;Prepare input files of Vl, V2, or V3;Analyze from files EEGV1, EEGV2, or EEGV3;Outfile file for plots */9.8FOURIER ANALYSIS OF AN ELECTROENCEPHALOGRAMFILE *fin,*fout;FILE *fopen();double twopi,yin,sigma,arg;double Vreal[MAX],Vimag[MAX],VFiltr[MAX],VFilti[MAX];/* Fortran compatibility; arrays are used from [l] */intnu,N,MAX1,j,it,choice,k,filter,maxf,maxk;char yn,yesno,wa,waFilter;void FFT();twopi = 8*atan(l.O);nu = 5; /* so */ N = 32; /* Number of time steps */printf("EEG FFT input & analysis:\n""Prepare input file? (y or n) :\n");scanf ("%s",&yn) ;if ( yn == 'y' )printf("Input EEG data set (1,2,3): ");scanf("%i",&choice) ;switch (choice)case 1: foutcase 2: foutcase 3: foutdefault:printf("\n !!}printf("\nInputfor ( j = 1; j= fopen ("EEGV1","w"); break;= fopen ("EEGV2","w"); break;= fopen ("EEGV3","w"); break;No such input file\n\n");exit(l);EEG data:\n") ;<= N; j++ )printf("\n%i: ", j);scanf("%lf",&yin) ;Vreal[j] = yin;fprintf(fout,"%lf\n",yin);fclose(fout); rewind(fout); /* Ready for reuse */printf("\nAnalyze which EEG data set (1,2,3)? ");scanf("%i",&choice) ;printf("\nWrite over (w) or Add on (a): ");scanf("%s",&wa):switch (choice) /* Open input & output files */case 1: fin = fopen("EEGV1","r");fout = fopen ("EEGout1",&wa); break;case 2: fin = fopen("EEGV2","r");369370DISCRETE FOURIER TRANSFORMS AND FOURIER SERIESfout = fopen ("EEGout2",&wa); break;case 3: fin = fopen("EEGV3","r");fout = fopen ("EEGout3",&wa); break;default:printf("\n !! No such input file\n\n");exit(l);for ( j = 1; j <= N; j++ ) /* Data input */fscanf(fin, "%lf\n",&yin) ;Vreal[j] = yin; Vimag[j] = 0;}it = 1; /* FFT from V back into V *//* Original V is destroyed */FFT(nu,N,it,twopi,Vreal,Vimag);printf("\nOutput of FFT; Vreal, Vimag\n");for ( k = 1; k <= N; k++ ) /* V has c-sub-k */fprintf(fout, "%lf %lf\n",Vreal[k],Vimag[k]);printf("\n%6.2lf%8.2lf",Vreal[k],Vimag[k]);/* Remaking EEG with filter options */yesno = 'y';tile ( yesno == 'y' )printf("\n\nFilter EEG; (y or n)? ");scanf ("%s" ,&yn) ; yesno = yn;if ( yesno == 'y' )printf("Choose filter type:\n");printf("1, no filter\n");2, Lanczos filter\n");printf("3, truncated filter\n");printf("scanf("%i",&filter);if ( abs(filter-2) > 1 )printf("\n !! No such filter\n\n");exit(l);printf("Write (w) or Add (a) FilterEEG%i\n',choice);scanf("%s",&waFilter) ;switch (choice) /* Gpen file for filtered EEG */case l:fout=fopen("FilterEEG1",&waFilter);break;case2:fout=fopen("FilterEEG2",&waFilter);break;case3:fout=fopen("FilterEEG3",&waFilter);break;9.8FOURIER ANALYSIS OF AN ELECTROENCEPHALOGRAMif ( filter == 3 ) /* Truncated filter */printf("Maximum frequency (0 to %i)? ",MAX-2);scanf("%i",&maxf); maxk = maxf+l;if ( maxk > MAX-1 ) maxk = MAX-l;for ( k = 1; k <= N; k++ )switch (filter)case 1: sigma = 1; break;/* No filter */case 2: /* Lanczos filter */if ( k == 1 ) { sigma = 1; break; >elsearg = twopi*(k-1)/(4*N);sigma = sin(arg)/arg; break;case 3: /* Truncated filter */if (k <= maxk) sigma = l;else sigma = 0; break;VFiltr[k] = sigma*Vreal[k];VFilti[k] = 0;it = -1; /* Remake voltages using FFT */FFT(nu,N,it,twopi,VFiltr,VFilti);printf("\nFiltered FFT; VFiltr,VFilti\n");for ( j = 1; j <= N; j++ ){fprintf(fout, "%lf %lf\n" ,VFiltr[j],VFilti[j]);printf("\n%6.2lf %8.2lf",VFiltr[j],VFilti[j]);}/* ends if yesno=='y' */}/* ends while yesno=='y' */printf("\nGoodbye from this brain wave");}/* ends main */FFT(nu,N,it,twopi,Vreal,Vimag) /* See Program 3, Section 9.7 */bitrev(j,nu) /* See Program 3, Section 9.7 */371372DISCRETE FOURIER TRANSFORMS AND FOURIER SERIESExercise 9.37(a) Code program EEG FFT Analysis, borrowing the functions FFT and bitrev from Program 9.3.
They are not listed in Program 9.4. Test the variousoptions to verify that the whole program is working correctly. For debugging, itis simplest to use a very small data set, nu = 2, so that the number of dataN = 4. This will minimize output and will obviate much retyping if you accidentally destroy your EEG data files.(b) Make another test of the FFT by using constant data, so that the Fourier coefficients for k > 0 should be zero. Also check that the inverse transform recovers the original data. nFrequency spectrum analysis of the EEGSince the EEG data in Table 9.1 are discrete and we are using the Fast Fouriertransform (FFT) method of analyzing the data, we will obtain from the program thediscrete Fourier transform (the DFT, Section 9.2) of the data from the time domainto the frequency domain.
Therefore, in our formalism in Section 9.2 the variable xis to be interpreted as the time variable t (in seconds), so that the variable k is interpreted as the frequency in hertz (Hz), that is, per second.The appropriate formulas for the DFT are those for the amplitudes ck, (9.6).The square of the modulus of each ck is proportional to the power contained in thatfrequency, as discussed under (9.8) for Parseval’s theorem.
Therefore, we definefor each frequency, k, the power, Pk, by(9.95)The program outputs the Fourier amplitudes ck, but you may easily also adapt it toproduce the power. In Figure 9.20 you can see the power spectrum from the FFTanalysis of data set V1 in Table 9.1.Now that you have seen a sample of the EEG analysis results, try your ownbrain waves on the data in Table 9.1.Exercise 9.38(a) Select an EEG-analysis output file EEGoutchoice, where choice = 1,2,or 3, and use it as input to a graphics application program in order to display thepower spectrum as a function of k, similarly to Figure 9.20.(b) Is the pattern of Fourier amplitudes that you obtained as a function of frequency, k, consistent with that of a normal adult human, as discussed in the introduction to the EEG? If you run more than one data set, compare them anddiscuss for what frequency range the analysis is relatively insensitive to noise.Recall that the average noise level increases from V1 to V2 to V3.
nNow that we have an understanding of the basic Fourier analysis of data, it istime to introduce some data massaging techniques.9.8FOURIER ANALYSIS OF AN ELECTROENCEPHALOGRAM373FIGURE 9.20 Power spectrum, Pk , as a function of frequency. k, for the EEG voltage inTable 9.1 having the lowest noise, namely V1(t).Filtering the EEG data: The Lanczos filterHere we explore how the reconstructed FFT can be modified by selectively downweighting various frequencies, usually the highest ones, that are usually associatedwith noise rather than with the desired signal.We use the code options for reconstructing V(t) by the inverse FFT.
Thus, inthe coding you see that the parameter passed to function FFT is it = -1, ratherthan it = 1, as used for the frequency analysis. It’s as simple as that for the discrete Fourier transform. Truncation-type filtering is a simple choice.Exercise 9.39This exercise consists of using the program segment in EEG FFT Analysis toreconstruct the voltage from the Fourier amplitudes computed in Exercise 9.38,but using a choice of filtering options, either none or truncation.(a) Use the filter type 1 (no filtering) to verify that the program reproduces theoriginal data, within the roundoff error of your computer and its arithmetic routines for cosines and sines.(b) Next, choose a maximum frequency, maxf , in EEG FFT Analysis. Thisfrequency may not be larger than the number of data minus one, but the programwill force this requirement. Run the program, then graph the output to compareit with the input data.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
