Introduction to Digital Signal Processing (DSP)

1. Introduction to Digital Signal Processing(DSP) Recent developments in digital computers open the way to this subject. The general block diagram of a DSP system is shownbelow: Bandlimit (pre- sampling)filter (antialiasing) A/D converter computer DSP Processor D/A converter Smoothing filter x(t) fs x(nTs) y(nTs) y(t) The input x(t) is an analogue signal (speech,video,…). 1-pre sampling filter: the signal is first bandlimited (for antialiasing*) using a LPF having cutoff fmax, where fmax is the maximum frequency content in x(t) ( to ensure that the highest frequency is within the bounds for which the signal can berecovered) 2-.A/D converter :where the bandlimited analogue signal is convertedinto digital x(nTs) using an analogue to digital converter (A/D) that consistsof: a sampler converts the continuos time signal into discrete time signal(with sampling frequency fs2fmax according to Nyquist theorem, whereTs=1/fs) a quantizer : assigning a binary value to the analogesamples. DSP: The digital signal (discrete) x(nTs) or simply written as x(n) is applied to a digital computer with a suitable interface card( sound card, video card,…). Now, inside this computer a program ( high level or low level) is written to perform any sort of signal analysis to x(n), like, linear or nonlinear amplification, convolution, correlation, digital filtering (LFP, HPF, BPF, BSF) and Fourier spectral analysis. The result of this processing is the digital signal y(nTs), or simplyy(n). D/A: This signal is then converted into analogue using a digital to analogue converter(D/A). 1

2. 5- smoothing filter: a smoothing filter is used to remove the staircase shape of y(n) to give the continuous output signaly(t). The software implementation of the signal analysis on x(n) giveshigh flexibility in design and processing. With the recent advances in digital technologies it is possible to implement very complicated systems on a single general purpose silicon chip. Classification of DSPsystems: According to the relation between the output y(n) and the input x(n), any DSP system can be classified according to thefollowing: 1-Linearity: A DSP system is called linear if the superposition theorem applies. Superposition means that sum of the input= sum of the output. Figure below illustrated that the system output due to the weighted sum inputs x1(n)+ x2(n) is equal to the same weighted sum of the individual outputs obtained from their corresponding inputs, that is,y(n)= y1(n)+y2(n). System x1(n) y1(n) System y2(n) x2(n) x1(n)+x2(n) System y1(n)+y2(n) Digital linearsystem Example1 if y(n)= 2x(n) which is a simple DSP system representing an amplifier withgain=2. This DSP system is linear since superposition applies here, if x(n)=x1(n)+x2(n),then: y(n)=2[x1(n)+x2(n)]=2 x1(n)+2x2(n)=y1(n)+y2(n), where y1(n) and y2(n)are the outputs due to x1(n) and x2(n) if they are applied one at atime. Example2y(n)=3e0.2x(n) which is a DSP system representing an exponential amplifier, now if x(n)=x1(n)+x2(n),then: y(n)=3e0.2[x (n)+x (n)]= 3e0.2x1(n) e0.2x2(n)≠ y (n)+y (n), this system isnonlinear. 1 2 1 2 Note : examples for the nonlinear systems are: y(n)= ln(x(n)),y(n)=sin(x(n)), y(n)=1/ x(n) , y(n)=(x(n))2 2

3. Causality: A DSP system is said to be causal if the present value of y [ i.e. y(n)] is not a function of some future values of x [i.e. x(n+i), i, +ve integer] or some future values of y [i.e. y(n+i), i, +ve integer]. The causal DSP system that is physically realizable is that system where y(n) is a function of present value of x(n) or some previous values of x [i.e x(n-i) ,i, +ve integer] or some previous values of y [i.e. y(n-i), i +veinteger]. In other words, a causal system is one where response does not begin before the input signal isapplied. Example1: if y(n)=3x(n)-2x(n-1)+y(n-3)+ex(n) which is not linear due tothe exponential term but it is causal since y(n) is not a function of some future terms of x ory. Example2: y(n)=x(n+1)-x(n)+3y(n+3) is not causal and is notphysically realizable since y(n) depends on some future values of xandy, [theseare x(n+1)) andy(n+3)]. Note : examples of the noncausal systems are: y(n)=x(n+1), y(n)=x(n2), y(n)=x(-n),y(n)=x│n│. Stability: A DSP system is said to be stable if the O/P is bounded for bounded I/P (bounded input gives bounded outputBIBO). Ex1: if y(n)=2x(n)-0.5x(n-1), and |x(n)|<G where G is finite, then|y(n)|<2G- 0.5G or |y(n)|<1.5G. hence if x(n) is bounded by G then y(n) is alsobounded i.e. the system isstable. Ex2: y(n)=ex(n-1) which is a nonlinear causal system, then if |x(n)|<G,then |y(n)|<eG and if G is finite, then y(n) is bounded, i.e., the system is stable. Stability of linearDSP system can be studied usually in terms of transfer function in z-domain, H(z)=Y(z)/X(z) Where the z-transforms of x(n) and y(n)are: X (z) x(n)z n n Y(z)y(n)zn n   If H(z) has pole(s) outside the unit circle (|z|=1), then the system is unstable, otherwise, it is a stable system. Poles on the unit circle gives criticallystable. 3

4. Imz unstable unstable Rez Stable 1 unstable unstable Ex: If y(n)=2x(n-1)—3x(n). Check if this system is stable or not. Solution: Taking Z-transform of both sides,get: Y(z)=2z-1X(z)-3X(z) (using shift property ofz-transform) OrY(z)/X(z)=H(z)=2z-1-3or: H(z)=(2-3z)/z which has a single pole at z=0(inside the unit circle), hence this system isstable. Ex: check the stability of the system:y(n)=x(n)+2y(n-1), Solution: As before, taking the z-transform of both sides and using the shift property, then: Y(z)=X(z)+2z-1Y(z), or Y(z)[1-2z-1]=X(z),then: H(z)=Y(z)/X(z)=1/(1-2z-1)=z/(z-2) which has a pole at z=2(outside the unit circle, hence this system isunstable. 4-Time Variant-Time Invariant DSPsystem: A time variant system is that system with time varying characteristics that depends on the time index n. For example, y(n)=n x(n) is a time variant system. In fact, this system is an amplifier having variable gain withtime. Take another example say y(n)=2e-x(n-1) which is a nonlinear, causal, stable system. This system is time invariant since its characteristics do not change with time indexn. Ex: Classify the following DSP system for linearity, causality, stability and timeinvariant: y(n)=2ex(n)-n x(n-2) +y(n-2). Solution: This systemis: -nonlinear due to the exponential term. -casual since y(n) does not depend on some future terms of x ory. 4

5. -unstable due to the term n x(n-2) which is unbounded with time index n even if x isbounded. -time variant due to the n x(n-2)term. Homeworks: classify the following DSP systems for linearity, causality, time-invariance andstability: a) y[n] = 3 x(n) – 4x(n-1) b) y[n] = 2 y(n-1) + x(n+2) y[n] = nx(n) y[n] = cos (x(n)) y[n] = log10(x(n)) y[n] =x[n]4 Input/output relations of the linearsystems: a) Analogue(continuous)systems(review): If h(t) is the impulse response, where t is the timeindex. H(w) is the transfer function which is the Fourier transform ofh(t). y(t)=x(t)  h(t) Where  is theconvolution: y(t) h(t) H(w) x(t) Y(w) X(w) t t y(t)x() h(t)dh()x(t)d 0 Also: Y(w)=X(w) H(w)and: |Y(w)|2 = |X(w)|2 |H(w)|2or 0 Gy(w)=Gx(w) |H(w)|2 where Gx(w) and Gy(w) are spectral densities of x and yrespectively.  Also output power y 2 (t) = | H (w) |2 G (w)dw x  5

6. b) Discrete (digital)systems: Here h(n) is the impulse response of the System,where n is the time index, and H(z) is its transferfunction: H(z)=Y(z)/X(z) And: y(n)=x(n)  h(n)or: y(n)x(n)h(nk)h(k)x(nk) h(n) H(z) x(n) y(n) X(z) Y(z) k k which is called the discreteconvolution. DSP systems are classified according to their impulse responses h(n)into: a) FIR (Finite Impulse Response): where h(n) has finite number of elements such as: h(n)={1,2,4,3,0,1} where the cursor indicates the position where n=0. h(n) 4 3 2 1 1 n b)IIR(Infinite Impulse Response): where h(n) has infinite number of elements such as: h(n)=(½)n u(n) where u(n) is the unit stepfunction: u(n)= 1 forn=0,1,2,3,… 0 elsewhere u(n) 1 n h(n) 1 0.5 0.25 0.125 n 6

7. Discrete convolutionMethods: 1- Graphical method: This includes the basic convolution steps :- reversing in time using k as time index, -shifting by n samples, -multiplication ofthe corresponding samples, -addition . x(k) 3 2 1 Ex: h(n)={1,-1,2}, x(n)={2,1,-1,3}. k -1 y(n)x(k)h(nk) k y(0)x(k)h(0k)=2+(-1)(-1)+(1)(3)=6 k h(-k) 2 1 y(1)x(k)h(1k)=2(-1)+(-1)(3)=-5 k k y(1)x(k)h(1k)=(2)(2)+(-1)(1)+(1)(-1)=2 k And so on, shifting of h(n) (left and right) until the overlapping between x(k) and h(n-k) disappears giving 0's at the outputy(n). -1 h(1-k) 2 1 k -1 h(-1-k) 2 1 k -1 7

8. 2-Tabular Method: This is a very simple method used for FIR systems with finite number of samples x(n). A rectangular table with N1 rows (no of elements in h(n)) and N2 columns(no of elements in x(n)), or visa versa, is arranged. Then, the cross multiplications are carried out. The sum of multiplications diagonally will give the values of y(n). Ex: Repeat previous example using tabular method. h(n)={1,-1,2},x(n)={2,1,-1,3}. 2 1 -1 3 1 -1 2 Theny(n)={2,-1,2,6,-5,6} Note that N=N1+N2-1=no of elements iny(n) =3+4-1=6 If O1 and O2 are positions of the cursors in h(n) and x(n)(from the left), then O=O1+O2-1=position of the cursor(position of y(0)) iny(n) =2+3-1=4 3- Add-overlapmethod: This is a modified method from the tabular method, when either h(n) or x(n) has large no of elements, then this can be divided into subsegments of smaller length. This helps to savememory. Ex: Find the discrete convolutionbetween: h(n)={1,-1,2} and x(n)={1,2,-1,3,4,-1,0,3} solution: Here x(n) is divided into segments of length 3 (say), with the last segment will be of length 2 (no problem). Hence previous tabular method is repeated 3 times: 8

9. 1 -1 2 1 2 -1 y1(n)={1,1,-1,5,-2} 1 -1 2 3 4 -1 y2(n)={3,1,1,9,-2} 1 -1 2 0 3 y3(n)={0,3,-3,6} Then, add y1, y2, and y3 , with y2 shifted to the left by 3 elements(length of each segment) and y3 shifted to the left by 6 elements,i.e.: y1 y2 y3 1 1 -1 5 -2 3 1 1 9 -2 0 3 -36 sum y(n)={ 1, 1, -1,8 , -1, 1, 9, 1, -3,6} with the cursor for y(n) at O=2+4-1=5 4- Matrixmethod: Here, [Y]=[A][h], where [Y]T=[y(n)]=output, [h]T=[h(n)] and the matrix [A] has N=N1+N2-1 rows and N2 columns(no of elements of h(n)). The 1st rowin [A] is the 1st element in x(n)(from the left) and the remaining elements are0's. The 2nd row in [A] is the 2nd element in x(n), then the 1st element in x(n) and the remaining elements are 0's. and so on until the last element in x(n) is entered at the N1th row( N1 is the no of elements in x(n)). After that 0's are applied in stead of the elements of x(n) until the last row atN=N1+N2-1. 9

10. Ex:h(n)={1,-1,2}, x(n)={2,1,-1,3} Solution: 2 2 0 0 1 2 0  1 1 2 1 3 0 2 0 0 1 2   1 1 2  1 1 3 1  0 0 3 1  2 1  0 1  [h]1 2  1 [ A] 3 1  6 then:[Y]= 3 =     2 5   0 0 1  6  0  3    Or y(n)={2,-1,2,6,-5,6} whereO=2+2-1=3 5-The Z-transform method: This is a general method used when either or both of h(n) and x(n) has infinite elements.(Note if both h(n) and x(n) has finite no of elements then use the tabular method since it is very simple). The procedure here is to take the z-transform of x(n) and h(n), then multiply to find Y(z) from which y(n) is found using inversez-transform. Ex: h(n)={1,-1,2}, x(n)=(½)nu(n). Solution: Note, since x(n) has infinite no of elements, then, we must use z- transformmethod. Now, taking the z-transforms of h(n) and x(n),then: H (z) h(n) z n =1. z1+(-1).z0+2.z-1=z-1+2z-1. n1 1  X(z)x(n)zn(0.5)nzn(0.5z1)n nn0 n0 or: (recall that in general Z[anu(n)]=z/(z-a)) X(z)1 z 10.5z1 z 0.5 ThenY(z)=H(z).X(z)=[z-1+2z-1][z/(z-0.5)] z.zz 2   z0.5 z0.5 z 0.5 And y(n)=(0.5)n+1u(n+1) –(0.5)n u(n) +2 (0.5)n-1u(n-1) Note: recall the shift property Z[x(n-a) u(n-a)]=X(z) z-a. 11

11. Ex: Find 1-(n-a)(n-b) 2-(n-a)u(n-b) 3- u(n-a)  u(n-b) using the z-transform and for real constants a & b. ((n) is the delta diracfunction). Solution: 1- Z[(n-a)]=z-a and Z[(n-b)]=z-b,then: Z-1[z-a z-b]= Z-1[z-(a+b)]=(n-{a+b}) Hence: (n-a) (n-b) =(n-{a+b}). 2- Z[(n-a)]=z-a and Z[u(n-b)]=z-b z/(z-1),then: Z-1[z-a z-b z/(z-1)]= Z-1[z-(a+b) z/(z-1)]= u(n-{a+b}) Hence: (n-a)  u(n-b) =u(n-{a+b}). 3- Z[u(n-a)]=z-a z/(z-1), Z[u(n-b)]=z-bz/(z-1), Z-1[z-a z/(z-1) z-b z/(z-b)]= Z-1[z-(a+b-1)z/(z-1)2] = (n-{a+b-1})u(n-{a+b-1}) which is a ramp function. (Note: recall thatZ[nu(n)]=z/(z-1)2 ). Hence: u(n-a)  u(n-b) =(n-{a+b-1})u(n-{a+b-1}) Homeworks: 1- find Z [sin0n]and Z[cos0n]. 2- find an u(n) u(n). 6-Frequency Response of the DSPsystem: If x(n) = Ae jo n (sampled sinusoid),then: y(n) h(k)x(n k) or: k h(n) y(n)h(k)Aejo(nk)Aejonh(k)ejok x(n) y(n) k k h(k)e  • jok H () If o which is the frequency response of h(n) at=o k of the input sinusoid. Or ingeneral, H()h(k)ejk |H()| () with magnitude and phasevalues. k So,the output y(n)can be easilyobtainedas: y(n) = x(n) |H()|() 11

12. Ex: For the DSP system shown, if x(t)=10cos(300t), find x(n) andy(n). A/D h(n)=(0.5)nu(n) D/A y(t) x(t) fs=1KHz Solution: fs=1000 Hz, then Ts=1/fs=0.001, t=nTs, then: x(nTs) = x(n) =10 cos(300(0.001n))=10cos(0.3n). hence, the input is a discrete sampled sinusoid with o=0.3 rad. To find y(n), we use the frequency response method since its much easier than the z-transform method. Then: y(n) = x(n) |H()|()  H() h(k)ejk(0.5)k ejk(0.5ej)k k k0 k0 Where 1 1  (recall geometric series r k0 1 ) 1r k Then: H (w) 10.5ejw Simplify: H () 10.5cos0.5jsin 1 1 |H()|  1.25cos 1cos0.25cos20.25sin2 ()  tan1 0.5sin 10.5cos For x(n)=10 cos(0.3n), then =o=0.3rad. 1 1.228and |H(0.3)| 1.25cos(0.3) (0.3) tan1 0.5sin(0.3) rad 10.5cos(0.3) 6 Then y(n) = 10 * 1.228 cos(0.3n) -/6=12.28cos{0.3n-(/6)}

13. DiscreteDeconvolution: To find h(n) ifbothx(n)andy(n)aregiven.x(n) h(n) y(n) Here, we usually use the z-transform method since it is valid for both the FIR and IIR DSPsystems. h(n)= Z-1[H(z)]=Z-1[Y(z)/X(z)] Ex: find h(n) if x(n)={ 1,-1,3} and y(n)={2,3,-1,17,-6} Solution: X(z)=1-z-1+3z-2 and Y(z)= 2+ 3z-1 -z-2 +17z-3- 6z-4. Then, using the long division: 2+5z-1 -2z-2 1-z-1 + 3z-2 2+ 3 z-1 – z-2 +17z-3-6z-4 2 - 2z-1 +6z-2 5z-1 – 7z-2 +17z-3-6z-4 5z-1 - 5z-2+15z-3 -2z-2 +2z-3-6z-4 -2z-2 +2z-3 -6z-4 0 0 0 Hence: H(z)= 2+5z-1 -2 z-2 or then, the system isFIR. h(n)={2,5,-2} and since there is noremainder, Ex: Find h(n) if x(n)=u(n) and y(n)=2u(n)-(0.5)nu(n). Solution: X(z)=z/(z-1) and Y(z)= [2z/(z-1)] –[z/(z-0.5)] 2zz  Then H (z) z1 z 0.5 2 z 1 2z 1z1z z z0.5 z0.5 z 0.5 z1 Hence h(n)=(0.5)nu(n).

14. DSP system Implementations a) FIRsystems: Here h(n) has finite no of elements: h(n)={h(0), h(1), h(2),……………..,h(m)} with (m+1)elements. H(z) = h(0) + h(1) z-1 + h(2) z-2 +…………….+ h(m) z-m and if : H(z)=Y(z)/X(z), then: y(n) = h(0) x(0) + h(1) x(n-1) + h(2) x(n-2) +…………….+ h(m)x(n-m) i.e., y(n) is obtained from x(n) by the weighted sum (weighted by h(n)) of the delayed samples of x(n). These delayed samples of x(n) are obtained using a tapped delay line with m-taps and with Ts time delay pertap. z-1Ts x(n) h(0) h(m) h(1)  y(n) Note that the FIR system is always stable, (no poles outside the unit circle), or we say that the FIR system is an open loop without a feedback and that is why it is called nonrecursive. b)-IIRsystem: For IIR system having m-zeros and r-poles,then: Y(z) a a z 1 a z 2 .......... . az m o 1 2 m H (z) X(z) 1bz1bz2........bzr 1 2 r Where, we can always set the 1st term at the denominator to unity by dividing with a suitable constant (say bo). Fromwhich: y(n)=ao x(n)+a1 x(n-1)+a2 x(n-2)+ ….+am x(n-m)- b1 y(n-1)- b2 y(n-2)-… ….- bry(n-r) Note that y(n) depends on present and previous values of x and on previous values of y.Hence: IIR system contains a feedback from the output to theinput. possibility of instability if some poles of H(z) lies outside the unitcircle.

15. 3- To implement the IIR system, then 2 tapped delay lines are required, one with m-taps for the x input and the other with r-taps for the feedback from the youtput.  y(n) x(n) b1 br Ex: Implement the DSPsystem: 42z2 H (z) 23z1z2z3 Solution: This is an IIR system. Dividing by 2 to set the 1st term at the denominator to unity,then: H(z)This needs a tapped delay line with 2 taps forx 11.5z10.5z20.5z3 2z2 and a tapped delay line with 3 taps fory.  x(n) y(n) z-1 z-1 z-1 z-1 z-1

16. Spectral Analysis of DiscreteSignals: a-periodic discretesignals: If x(n) is a periodic signal having a period of N samples, i.e. x(n) is repeated every N samples, then the frequency content of x(n) (spectral analysis) is obtained using discrete Fourier series (DFS),where: • j2nk • N 1 N N1 x(n)e X(k) N- 1k0 n0 where n is the time index and k is the frequency index, i.e.: X(0)=dccomponent. X(1)=1st harmoniccomponent. X(2)=2nd harmonic component and so on up to the N-1 harmonic component X(N-1). j 2nk N1 x(n)X(k)e N The inverse DFS is givenby: N- 1n0 k0 Ex: Find the DFS of the signals x(n)shown: Solution: x(n) n Note that x(n) is periodic with period N=4 samples and is defined for full period by: x(n)={2,-1,3,-3}.Then: j 2n(0) 4 x(0)x(1)x(2)x(3)21331 4 1 3 X(0)x(n)e   =dcvalue 4 4 4n0

17. j 2n(1) 4 j3 2] 1 1 3 X(1)x(n)e • j • 2 • j [x(0) x(1)e x(2)e x(3)e 4 4n0 X(1)=0.25[2-1(-j)+3(-1)-3(j)] = 0.25(-1-2j) =-0.25-0.5j Similarly, we find X(2)=9/4 and X(3)=-0.25+0.5j Notethat X(3)=X*(1) (*conjugate) It should be noted, however, that the discrete behavior of DFS directly gives the discrete line spectrum of Fourier series analysis with k index being the harmonic number. Also this frequency index is for N- 1k 0 , i.e. for N- sampled periodic signal, the harmonic analysis is done up to(N-1). b- aperiodic or random discretesignals: If x(n) is aperiodic or random, such as pulses, speech, video,…….. when sampled at fs then spectral analysis is done for finite segment (length) of such signals, i.e. N samples are taken from such a signal ( this corresponds to NTs time interval from x(t)). Here, we use the discrete Fourier transform(DFT). x(t) x(N-1) x(0) N samples If x(n)={x(0), x(1),x(2),………….x(N-1)} is the sampled signal from x(t), then the DFT of x(n) is givenby: • j2nk • N 1 N x(n)e n0 N1 X(k) N- 1k0 j 2nk N1 Andthe inverseDFT(IDFT)is:x(n)X(k)e N k0 Theoretically, x(t) has continuous spectrum, but due to sampling at fs, then this will give a discrete approximation to this continuous spectrum byX(k).

18. The relation between frequency in Hz and the frequency index k is given by the frequency resolution between successive values of k, this is givenby: frequncyresolutionfs N Such that the discrete behavior of |X(k)| approaches the continuousspectrum asN. |X(k)| k N/2 N N-1 fs/2 fs fHz fs/N Note also that there is a conjugate and even symmetry around (N/2), and that the spectrum repeats itself after k>N. Ex: Speech signal is sampled at 8000Hz. Find the min number of samples to the DFT analysis such that frequency analysis is done with a resolution of 20Hz. Solution: for 20 Hz resolution, then 20=fs/N=8000/N, then N=400 samples, i.e.: x(n)={x(0),x(1),x(2),………….x(399)} and: • j2nk • 400 1 399 x(n)e X(k) 399  k0 400n0 Note that above summation needs 400*400=160000 complex multiplications to find the whole spectrum (399  k 0), or in general DFT requires N2 complex multiplications. For such large number of calculations, we usually use the digitalcomputer.

19. Ex: Find the DFT of the sampled sequence representing a sampled pulse: x(n)={1,1,1,0,0,0}. x(n) 1 1 1 1 2 n X(0)=(1/6)(1+1+1)=0.5=dcvalue |X(0)|=0.5 j 2n(1) 6 1 5 X(1)x(n)e =(1/6)[1+e-j/3+e-j2/3]=(1/6)[1+0.5-j0.866+(-0.5)-j0.866 X(1)=(1/6)[1-j1.732] ,|X(1)|=0.333 6n 0 j 2n(2) 6 1 5 X(2)x(n)e =(1/6)[1+e-j2/3+e-j4/3]=(1/6)[1-0.5-j0.866-0.5+j0.866] 6n0 X(2)=0,|X(2)|=0 j 2n (3) 1 5 X(3)x(n)e 6 =(1/6)[ 1+e-j+e-j2]=(1/6)[1-1+1] 6n0 X(3)=1/6,|X(3)|=0.1666 j 2n(4) 1 5 X(4)x(n)e 6 =(1/6)[1+e-j4/3+e-j8/3]=(1/6)[1-0.5+j0.866-0.5-j0.866] 6n0 X(4)=0,|X(4)|=0 j 2n(5) 1 5 X(5)x(n)e 6 =(1/6)[1+e-j5/3+e-j10/3]=(1/6)[1+0.5+j0.866-0.5+j0.866] 6n0 X(5)=(1/6)[1+j1.732],|X(5)|=0.333 0.5 1/3 1/3 1/6 1 2 3 4 5 Note the even symmetry and the complexconjugate Note also that the energy if computed from time domain,then:

20. Ex2(n)1113 n0 And if computed from frequency domain, then: E=N|X2(k)|6[0.25(1/9)(1/36)(1/9)]3E fromtimedomain k0 5 N1 N1 N1 x2(n)N|X2(k)| n0 k0 Or we say in generalthat: Fast Fourier Transform(FFT) This is a fast method to find the DFT. For FFT base-2 (radix-2), then N must be a power of 2, i.e. N=2r ( 4,8,16,32,64,128,256,512,……). If N≠2r, then 0's are added to complete the sequence to the nearest 2rvalue. Algorithm for FFT radix-2 decimation intime: • j2nk • N 1 N1 x(n)e 1)First, we know thatX (k )  , where the term (1/N) isa Nn0 • j 2nk N1 X(k)x(n)e N 2)scale factor that can be omitted,then n0 2/N If we define W =e-j2/N, to ease notation,then: N N1 N  neven N  nodd nk nk nk N X(k)x(n)W x(n)W x(n)W -j2/N e n0 for n even, then n=2r and for n odd, then n=2r+1, wherer=0,1,2,3,….. N1 N1 2 2 N  2rk (2r1)k X(k)x(2r)W x(2r1)W 3) N r0 r0 WN2=e-j4/N=e-j2/(N/2)=WN/2,then: Notethat N1 N1 2 2 N/2 NN/2 rk k rk X(k)x(2r)W W x(2r1)W r0 r0 X(k)=G(k)+WNk + H(k)………….(1) where G(k)=DFT of the even numbered samples (N/2 samples) H(k)=DFT of the odd numbered samples (N/2 samples). 21

21. Equation(1) states that the DFT of N samples can be obtained in terms of the DFT of the even and odd numbered N/2 samples. And if eq(1) again used to find the DFT's of the N/2 samples in terms of the DFT's of N/4 samples. This process is repeated many times ( exactly r=log2N times) until we end up with the DFT of 2 samples. DFT of 2-samples: • j 2nk 1 Ifx(n)={x(0),x(1)},then:X(k)x(n)e n0 X(0)=x(0)+x(1) X(1)=x(0)-x(1) Using signal flow representation,then: 2 fromwhich: x(0) X(0) x(1) X(1) Timedomain freqdomain where the unity path gain is notmarked. The 2-point (samples) DFT is called "Butterfly''. Note: for N=2r point FFT, then segmenting into even and odd is done for r times and the number of complex multiplications in FFT will be r.N or: Number of complex multiplications in FFT= Nlog2N. Note that, the great advantage of using FFT is when N islarger. 21

22. The process of finding say 64-point FFT is started with 2-point FFT(Butterfly). This is then used using eq(1) to find the 4-point FFT, and again, this 4-point FFT is used to find the 8-point FFT and soon. Ex: Draw the signal flow graph of the 4-point FFT, then use to it to find the spectrum of the sequencex(n)={1,-1,2,3}. Solution: X(k)=G(k)+WkH(k) 3  k0. 4 Theeven numbered samples are x(0) and x(2) The odd numbered samples are x(1) and x(3) X(0)=G(0)+W 0H(0) 4 X(1)=G(1)+W 1H(1) 4 Mirror X(2)=G(2)+W 2H(2) 4 since G(2)=G(0)and H(2)=H(0) where 2-point DFT has a period of 2samples,then: X(2)=G(0)+W 2H(0) 4 X(3)=G(3)+W 3H(3) 4 And again G(3)=G(1)and H(3)=H(1),then: X(3)=G(1)+W 3H(1) 4 G(0) 2-point DFT Butterfly x(0) X(0) x(0) x(2) x(1) x(3) 2-point DFT Butterfly G(1) x(2) X(1) 2-point DFT Butterfly x(1) X(2) 2-point DFT Butterfly H(0) x(3) X(3) H(1) G(0)=3 X(0)=5 x(0) X(1)=-1+j4 x(2) G(1)=-1 H(0)=2 X(2)=3-2=1 x(1) H(1)=-4 x(3) X(3)=-1-j4 22

23. Check: x2 (n) 1+4+1+9=15 =energy from time domain 1|X(k)|2(25+17+1+17)/4=15=energy fromfrequency domain N Ex: Draw the signal flow graph of 8-pointFFT. Solution: The 8=point FFT uses the 4-point FFT obtained in previous example. The sequence of sample decimation(splitting into even and odd numbered samples) is obtained using mirror image of 3-bit dataas: x(0) x(4) x(2) x(6) x(1) x(5) x(3) x(7) X(0) 4-point FFT X(1) X(2) X(3) X(4) 4-point FFT X(5) X(6) X(7) Time domain Frequency domain Mirror

24. X(0)=G(0)+W80H(0) x(0) X(1)=G(1)+W81H(1) x(4) X(2)=G(2)+W82H(2) x(2) X(3)=G(3)+W83H(3) x(6) X(4)=G(4)+W84H(4) x(1) x(5) X(5)=G(5)+W85H(5) X(6)=G(6)+W86H(6) x(3) X(7)=G(7)+W87H(7) x(7)