EE406 Discrete-Time Signal Processing
Fast Fourier Transform (FFT) Algorithms
Linear Convolution vs Circular Convolution in the DFT The linear convolution of two sequences x1 [ n] of N 1 -point and x 2 [ n] of N 2 -point is given by x3 [ n] = x1 [ n] ∗ x 2 [ n] =
N 1 −1
∞
∑ x [k ] x [n − k ] = ∑ x [k ] x [n − k ] 1
2
k = −∞
1
2
k = 0
x3 [ n] is a ( N 1 + N 2 − 1) -point sequence. If we choose N =max( N 1 , N 2 ) and compute a N -
point circular convolution x1 [n] ⊗ x 2 [n] , then we obtain a N -point sequence, which is obviously different from x3 [ n ] . If we choose N = ( N 1 + N 2 − 1) and perform a ( N 1 + N 2 − 1) -point circular convolution, the result becomes
N −1 x 4 [ n] = x1 [ n] ⊗ x 2 [n] = ∑ x1 [ m] x 2 ([n − m]) N R N [ n] m=0 ∞ N −1 = ∑ x1 [m] ∑ x 2 [n − m − rN ] R N [n] r = −∞ m =0 ∞ N −1 = ∑ ∑ x1 [m] x 2 [n − m − rN ] R N [n] r =−∞m =0 ∞ = ∑ x3 [n − rN ] R N [n] r =−∞ Thus the circular convolution is an aliased version of the llinear inear convolution. Since x3 [ n] is a N = ( N 1 + N 2 − 1) -point sequence, we have x 4 [ n] = x 3 [n]
0 ≤ n ≤ N − 1
which means that there is no aliasing in the time domain. If we make both x1 [ n] and x 2 [n] a N = ( N 1 + N 2 − 1) -point sequence by padding an appropriate number of zeros, then the circular convolution is identical to the linear convolution.
In order to use the DFT for linear convolution, we must choose N properly properl y. However, when N is large, there is an immense requirement on memory. If N is chosen to be less than the required value, an an error will be introduced. Suppose we choose choose N such that max( N 1 , N 2 ) ≤ N < ( N 1 + N 2 − 1) Since
∞ x 4 [ n] = ∑ x3 [ n − rN ] R N [ n] , r =−∞
the error introduced is then
ZA-041
Page - 8
EE406 Discrete-Time Signal Processing
Fast Fourier Transform (FFT) Algorithms
∞ e[ n] = x 4 [n ] − x3 [n ] = ∑ x3 [n − rN ] R N [n] . r ≠0 r =−∞ In general x1 [n] and x 2 [n] are causal, thus x3 [ n] is also causal. As max( N 1 , N 2 ) ≤ N , only two terms r = ±1 remain, and e[ n] = [ x3 [ n − N ] + x3 [ n + N ]] R N [ n] x3 [n − N ] = 0, 0 ≤ n ≤ N − 1 e[ n] = x3 [n + N ]
then
0 ≤ n ≤ N − 1 .
Thus even though we choose max( N 1 , N 2 ) ≤ N , the error sample at n is the same as the linear convolution N samples away. Now the linear convolution will be zero after ( N 1 + N 2 − 1) samples. This means that the first few samples of the circular convolution are in error, while the remaining ones are the correct linear convolution values. Example: Let x1 [ n] and x 2 [ n] be the two 4-point sequences given by { x1 [ n]} ={1, 2, 2, 1} and { x 2 [ n]} ={1, −1, 1, −1} a. Determine the linear convolution x3 [n] = x1[ n] ∗ x2 [n] . b. Compute the circular convolution x 4 [ n] = x1 [ n] ⊗ x2 [ n] for N = 6, 5, and 4, and verify the error relations in each case. [Solution] The linear convolution is x3 [n ] ={1, 1, −1, −2, −1, 1, 1}.
When N =6, 6-point sequence x 4 [ n] = x1 [ n] ⊗ x2 [ n] ={2, 1, −1, −2, −1, 1}
Therefore e[n] = {2, 1, −1, −2, −1, 1} − {1, 1, −1, −2, −1, 1} = {1, 0, 0, 0, 0, 0} = x3 [ n + 6]
0≤n≤5
When N =5, 5-point sequence x 4 [ n] = x1 [ n] ⊗ x2 [ n] ={2, 2, −1, −2, −1}
and ZA-041
Page - 9
EE406 Discrete-Time Signal Processing
Fast Fourier Transform (FFT) Algorithms
e[n] = {2, 2, −1, −2, −1} − {1, 1, −1, −2, −1} = {1, 1, 0, 0, 0} = x3 [ n + 5] .
0≤n≤4
Finally when N = 4, the 4-point sequence is x 4 [ n] = x1 [ n] ⊗ x2 [ n] ={0, 2, 0, −2} e[n] = {0, 2, 0, −2} − {1, 1, −1, −2} = x3 [n + 4]
0≤n≤3
Thus when N = max( N 1 , N 2 ) is chosen for circular convolution, the first ( M –1) samples are in error where M = min ( N 1 , N 2 ) . This result is very useful for the implementation of long convolutions in the form of block processing.
Block Convolution of the DFT for Spectral Analysis For an infinite-length sequence, a practical approach is to partition the infinite-length input sequence into smaller blocks, process each block using the DFT, and finally assemble the output sequence from the outputs of each block. This procedure is called a block convolution. Overlap-save method for block convolution
Consider an input sequence x[n] being partitioned into N -point sequences and the impulse response of the filter is a M -point sequence where M < N . If we simply partition x[n] into nonoverlapping blocks, then the resulting output sequence will have intervals of incorrect samples. To correct this problem, we can partition x[n] into overlapping blocks with the previous one by exactly ( M –1) samples, save the last ( N–M +1) output samples, and finally concatenate these outputs into a sequence. To correct for the first ( M –1) samples in the first output block, we set the first ( M –1) samples to zero. This procedure is called overlap-save method. Example: Let x[n]=(n+1), 0 ≤ n ≤ 9 and h[ n] = {1, 0, − 1} . Calculate the linear convolution with N =6 by using the overlap-save method. [Solution] Now M = 3, N = 6 and we will need ( M –1)=2 zeros at the beginning. We need to have 3 blocks with x1 [ n] = {0, 0, 1, 2, 3, 4} x 2 [ n] = {3, 4, 5, 6, 7, 8} x3 [n] = {7, 8, 9, 10, 0, 0}
ZA-041
Page - 10
EE406 Discrete-Time Signal Processing
Fast Fourier Transform (FFT) Algorithms
We pad x3 [ n] with two zeros since x[n] is running out of data. Now y1 [n] = x1 [ n] ⊗ h[ n] = {–3, –4, 1, 2, 2, 2} y 2 [n] = x 2 [ n] ⊗ h[ n] = {–4, –4, 2, 2, 2, 2} y3 [ n] = x3 [ n] ⊗ h[ n] = {7, 8, 2, 2, –9, –10}
y[n] = x[ n] ⊗ h[n] ={1, 2, 2, 2, 2, 2, 2, 2, 2, 2, –9, –10}.
Overlap-add method for block convolution
There is another method called overlap-add method for block convolution. In this method, the input sequence x[n] is partitioned into nonoverlapping blocks and convolved with the impulse response. Let x[n] be a long sequence of length ML with M , L>>0, and h[n] is a Lpoint sequence. First divide x[n] into M nonoverlapping blocks with each length L, and
x[n], kM ≤ n ≤ (k + 1) M − 1 x k [ n] = otherwise 0, Then y[n] = x[ n ] ∗ h[n ] =
M −1
M −1
∑ x [n] ∗ h[n] = ∑ y [n] . k
k =0
k
k =0
where y k [ n] is a (2 L – 1)-point sequence. Example: Recalculate the convolution in the previous example by using overlap-add method. [Solution] Now L = 3, x1 [ n] = {1, 2, 3} x 2 [ n] = {4, 5, 6} x3 [n] = {7, 8, 9}
x 4 [ n] = {10, 0, 0} y1 [n] = x1 [ n] ∗ h[n] = {1, 2, 2, –2, –3} y 2 [ n] = x 2 [ n] ∗ h[n] = {4, 5, 2, –5, –6} y3 [ n] = x3 [ n] ∗ h[ n] = {7, 8, 2, –8, –9}
y 4 [ n] = x 4 [ n] ∗ h[n] = {10, 0, –10, 0, 0} y[ n] = x[n ] ∗ h[n ] = {1, 2, 2, –2, –3} + {4, 5, 2, –5, –6} + {7, 8, 2, –8, –9} + {10, 0, –10, 0, 0} = {1, 2, 2, 2, 2, 2, 2, 2, 2, 2, –9, –10} ZA-041
Page - 11