Computing Fourier Series and Power Spectrum with MATLAB

1 ComputingFourierSeriesand PowerSpectrumwithMATLABBy IntroductionFourierseriesprovidesan alternateway of representingdata:insteadof represent-ing the signalamplitudeas a functionof time,we represent the signalby how muchinformationis containedat di erent you ever watchedthe blink-ing lights on a stereoequalizerthenyou have seenFourieranalysisat represent whetherthe musiccontainslots of bassor ,a French Mathematicianwhoonceserved as a scienti cadvisertoNapoleon,is creditedwiththe discovery of the resultsthatnow bear his important in dataacquisitionjustas it is in you might want to boost the Power of the basson your stereoyou might wantto lterout highfrequencynoisefromthe nearby radiotowers in Needhamwhenyou are conductinga lab experiment.

2 Fourieranalysisallows you to will describe someof the basicsof Fourierseriesand will showyou how you can it easyto translatea signalfromthe timedomainto the frequencydomain,one must understandhow to interpretthe datain the TheMathPartThissectionwill not providerigorousderivationsof Fourierseriesand it's proper-ties;thosetopicswill be coveredlaterin one of your the generalnatureof the Fourierseriesand how transformfunctionsfromthetimedomainto the Fourierseriestakes a signaland decomposesit into a sumof sinesand cosinesof di erent have a signalthatlast for 1 second,0<t <1, we conjecturethatcan represent thatsignalby the in niteseriesf(t) =a0+1Xn=1(ansin(2 nt) +bncos(2 nt))( )wheref(t) is the signalin the timedomain,andanandbnare unknown coe cientsof the ,n, has unitsof Hertz(Hz)=1=sandcorrespondsto thefrequencyof the wave.

3 We will not prove in this sectionthatthe seriesconvergesorthatthis seriescan providean accuraterepresentationof any will 'prove'by exampleto show thatthisseriesactuallyworksandcan will comelaterin your we de nedthe Fourierserieson a timeinterval 0< t <1, but we can alsodo a make the above equationworkTEX Paper2 FourierSeriesfor any timeinterval. We willdiscussdi erent timeintervals later,but willusethe one secondinterval for convenienceat thispoint. Alsonotethatthe Fourierrepresentationis formallyperiodic, the beginningof the cyclewill always equaltheend ( (t= 0) =f(t= 1)).To ndthe coe cientsanandbnyou can use the be derived in the appendix,but we will just statethe resulthere:Z10sin(2 nt)sin(2 mt)dt= 0; n and m are integers,and n6= m( )Z10sin(2 nt)sin(2 nt)dt= 1=2; n is an integer( )Z10cos(2 nt)cos(2 mt)dt= 0; n and m are integers,and n6= m( )Z10cos(2 nt)cos(2 nt)dt= 1=2; n is an integer( )Z10cos(2 nt)sin(2 mt)dt= 0; n and m are integers( )Thesepropertiesmake ndingthe coe cients quiteeasy.

4 You can multiplybothsidesof bysin(2 mt) andintegratethe expressionover the interval0< t <1. Integratingboth sidesof an equationis no di erent thanmultiplyingboth sidesby a constant: as longas you performthe sameoperationto both sidesof the equationthe equality will f(t) =a0+1Xn=0(ansin(2 nt) +bncos(2 nt))!sin(2 mt)dt:( )Applyingthe propertiesstatedabove, we know thatthe for all termsin the sumwherem does not equaln the resultof the integralwill be alsoknowthatany termsthatinvolve the productof sineand cosinewill be the sumwill be the sine termwheren=m. Therefore, will reducetoZ10f(t)sin(2 nt)dt=an2:( )Likewise,we can bycos(2 mt) and integrateover 0< t <1to obtainthe analogousresult:Z10f(t)cos(2 nt)dt=bn2:( )Thecoe cienta0is even easier:we onlyneedto over theinterval 0< t <1 to (t) =a0:( )Notethatwe have usedthe fact thatall the sine and cosinetermsintegrateto zeroover the domainof we have at our disposal , , & us to computethe unknown coe SomeexamplesTheeasiestexamplewouldbe to setf(t) =sin(2 t).

5 Withouteven performingthe calculation(simplyinspect )we know thatthe Fouriertransformshouldgive usa1= 1 and all othercoe cients shouldbe check thatthisworks,insertthe test functionf(t) =sin(2 t) into and to see (2 t)sin(2 nt)dt=an:( )Usingthe integralpropertiesgiven previously, it is clearto see thatall a's are zeroexceptfora1= 1: the resultwe thebncoe cients are zeroas wellbecausethe integralofsin(2 t)cos(2 nt) will always be is alsoeasyto check thatif f(t) is a constant thenall the coe cients exceptfora0will be zero, 'stry computinga Fourierseriesfor a squarewave signalthatis on for half thetimeinterval and the o for half, (t) = 1; 0< t <1=2 andf(t) = 0; 1=2< t < get the coe cients we and (2 nt)dt=an:( )Carryingout the integralover halfthe cycle(theintegralover the interval 1=2<t <1 will always be zerosincethe signalis zeroduringthis time) 2cos(2 nt)2 n 1=20=an( )and evaluatingatt= 0 andt= 1=2 yields1 cos( nt) n=an:( )Thereforewhenn is odd,an=2 n( )and whenn is evenan= 0( )Carryingout the sameoperationfor the coe (2 nt)2 n 1=20=bn( )and evaluatingatt= 0 andt= 1=2 yieldssin( nt) n=bn;( )which, for all n, is simplybn= 0( ) (t)n=20 n=2,000 Figure1.

6 Fourierreconstructionof a squarewave. We see a lot of ringingin the seriesuntil we includemany points into the also needthe coe cient,a0, which is obtainedby overthe timeintervala0=Z1=20dt= 1=2:( )Now we can plotthe FourierseriesrepresentationusingMATLABan d see howthe seriesdoes at reproducingthe we plotthe rst20 termsin thesum,we see the generalshape of the originalfunctionbut we see a lot of 'ringing'.Aswe plotmoreterms,we see the originalfunctionis represented quiteaccurately. IngeneralFourierseriescan reconstructa signalwitha smallnumber of modes if theoriginalsignalis highfrequencycomponentsto constructthe reference,the MATLAB code thatgeneratedthis gureis given below:N = 2000;x = [0:100]/100;f = ones(1,101)*1/2;fori = 1:2:Na = 2/pi/i;f = f+ a*sin(2*pi*i*x);endplot(x,f)4.

7 FourierTransformof DiscreteDataWhenwe are workingwithdatathatwe have acquiredwiththe DAQ cardwe havediscretedatapoints, not an analyticalfunctionthatwe can turnsout thattakinga Fouriertransformof discretedatais doneby simplytakinga discreteapproximationto the integrals( , , ). 1 (t)t1t2t3t4t5t6t7t8A1A2A3A4A5A6A7 Figure2. Schematicof trapezoidalruleappliedto the functionf(t) =sin(2 t) sampledat 8 points. Thetrapezoidalrulesimplycomputesthe areaof each block and be computedusingthe a functioninto several trapezoidsand sumsthe area,see Figure2In Figure2 we see thatthereare seven blocks (Aj) and 8 datapoints (tj). Tocomputethe area,Aj, of thejthtrapezoidalblock we useAj=f(tj) +f(tj+1)2(tj+1 tj)( )For now, let'sassumethatthe pointstjare equallyspacedsuch thattj+1 tj= t.

8 Thereforethe totalareaisA= t f(t1) +f(t2)2+f(t2) +f(t3)2+:::f(t6) +f(t6)2+f(t7) +f(t8)2 ;( )which generalizestoZf(t)dt= t0@f(t1)=2 +f(tn)=2 +N 1Xj=2f(tj)1A( )Thetrapezoidalruleprovidesa generalmechanismfor computethe Fouriercoe cient,an, you simplyneedto applythe trapezoidruleto the functionmultipliedby the functionsin(2 nt), speci cally,an= t[sin(2 nt1)f(t1) +sin(2 ntN)f(tN) + 2N 1Xj=2sin(2 ntj)f(tj)]( )ExerciseFS 1:Write a generalMATLAB functionthattakes a twovectors,x andy as input,assumesy is a functionof x andcomputes6 FourierSeriesthe integralRx(N)x(1)y(x)dxusingthe trapezoidalrule ( ).SomeMATLAB functionsthatmay be useful,thoughsurelynot necessary, arelength,sum,y(end),diff, andunique.

9 Assumesthatthepoints in x are equallyspaced- your functioncan make this assumptionas well. Test the numericalintegralon the functionssin(2 x) andx2on the interval of 0< x <1. Does the numericalapproximationmatchthe analyticalresult?Afteryou writeyour function,which you cannameintegral, you shouldbe ableto go to the MATLAB promptandissuethe followingcommands:x = [0:100]/100;y = sin(2*pi*x);to de nex and y. You shouldthenbe ableto typeintegral(y,x)withno semi-colonto get the valueof the integralreturnedat the the expressionfor the numericalintegralwherethe points in x are not functiontowork for a generalspacingin x. Test your functionusingan x de the approximationis de nedasthe di erencebetween the 'true'analyticalintegraland the errorin the integralforsin(2 x) in a rangeof 10to 10,000datapoints.

10 Ploton a log-logplotthe errorvs. the number 2:Write a generalMATLAB functionthattakes a discretefunctionf(tj) anda frequencynas inputandcomputesthe discreteFouriercoe cient Assumethatthe inputfunctionspansfrom0< t <1 andthatthe datapoints are equallyspacedintimeon that0< t <1 interval: t= 1=(N 1) where tis thespacingin timebetween points and N is the number of points usedtorepresent the functionf. Createa discretefunctionusing20 points thatcorrespondsto the squarewave thatwe usedin Figure1. For example:f= [ones(1,10)zeros(1,10)];de nesthe squarewave. Numericallycomputea few di erent Fouriercoe cients (bothanandbncoe cients)withyour MATLAB function:try n=0,1,2,3and 4. Show thatthe nu-mericalresultis in good agreement withthe analyticalexpressionforthe squarewave, your functionto computethe Fourierco-e cients for all frequencies,n.