Monday, 30 March 2015

fft - DFT as convolution question


I have tried to make this question as readable and consistent as possible. The short of it, is that I am trying to ascertain how one gets from the math equation shown, (which I understand), to the correct implementation. (Which I also understand). However, how/why one gets from one to the other seems unclear to me.


Here is my setup:


We have an $M$ length 'original' signal $x[m]$: (Equation 1)


$$ x[m] , \space m = 0,1,..M-1 $$


We wish to take a 'zoom-DFT' of it, where the DFT length will be, $N$, where $N > M$.


We thus now have an $N$ length zero-padded new signal, lets call it $x_p[n]$, where: (Equation 2) $$ x_p[n], \space n = 0,1,...N-1 $$


Our discrete time variable is now $n = 0,1,...N-1$, and our discrete frequency variable is $k = 0,1,...N-1$.



The $N$-length DFT of a signal is defined as such: (Equation 3)


$$ X[k] = \sum_{n=0}^{N-1} x_p[n]\space e^{-j\frac{2\pi nk}{N}} $$


If the $nk$ in the DFT equation is written as $\frac{-(k-n)^2 + n^2 + k^2}{2}$, then the DFT can now be written as: (Equation 4)


$$ X[k] = e^{-j\frac{\pi k^2}{N}} \space \sum_{n=0}^{N-1} x_p[n]\space e^{-j\frac{\pi n^2}{N}} \space e^{j\frac{\pi (k-n)^2}{N}} $$


This expresses the full DFT, (i.e, the DFT from 0 to nyquist), as a convolution. To evaluate the 'Zoom-DFT', where we want to interrogate only frequencies between $f_1$ and $f_2$, and so the following minor modification is made: If we let $F_{\Delta} = \frac{(f_2 - f_1)}{f_s}$, and $F_1 = \frac{f_1}{f_s}$, then the zoom-DFT, $X_{z}[k]$ becomes: (Equation 5)


$$ X_z[k] = e^{-j\frac{\pi F_{\Delta} k^2}{N}} \space \sum_{n=0}^{N-1} x_p[n]\space e^{-j\frac{\pi n^2}{N}} \space e^{-j 2\pi \space F_1 n}\space e^{j\frac{\pi F_{\Delta} (k-n)^2}{N}} $$


We can then simplify, if we let: (Equation 6)


$$ a[u] = x_p[u]\space e^{-j\frac{\pi u^2}{N}} \space e^{-j 2\pi \space F_1 u} \\ b[u] = e^{j\frac{\pi F_{\Delta} u^2}{N}} $$


, where $u$ is just some general dummy variable, (you will see why I picked it later). Anyway, now finally, we simply have: (Equation 7)


$$ X_z[k] = b^*[k] \space \sum_{n=0}^{N-1} a[n] \space b[k-n] $$



This, as we are told here, and here, page 183, is how the DFT can be re-written as a convolution.


Great! So far so good. The problem comes in interpretation vs implementation.


The problem:


My interpretation of (Equation 7), is that $a[u]$ is of length $N$, and $b[u]$ is also of length $N$, where we evaluate then for $u = 0,1,...N-1$. Then, we do a circular convolution of the series $a[u]$ and $b[u]$, to attain a series also of length $N$, (which is our desired DFT length), before we finally perform an element-by-element post-multiply with another series, also of length $N$. I dont so much care about post-multiply, that is not the point of this question. The point here is that $a[u]$ is length $N$, $b[u]$ is also length $N$, and $u$ for both of them is evaluated for $0 \leq u \leq N-1$, and then we do a circular convolution, modulo-$N$.


If you however now take this interpretation and implement it, you will not get the right answer. In fact, you will only ever get the right answer, if the $u$ variable for evaluation of $a[u]$ is taken from $u = 0,1,...N-1$, (same as above), BUT, the $u$ variable for evaluation of $b[u]$ is taken from $-(M-1) \leq u \leq N-1$. So now the series $b[u]$ is not length $N$ anymore, it is of length $N+M-1$. Then, if we do a modulo-($N+M-1$) circular convolution, and extract the indicies from $M$ to $M+N-1$, we will get the right answer.


...How, from equation-1, does one possibly get this 'correct' interpretation though? This, I cannot seem to explain, however I understand everything else.


FWIW, I have sat down and done the computation 'by hand', so I can 'see' why $b[u]$ MUST be taken from $-(M-1) \leq u \leq N-1$. (It would otherwise be impossible to compute the DFT for $k = 0,1,...N-1$). However I do not understand how one goes from (Equation 7), to this conclusion. Right now, to me, EQ-1 just looks like a garden variety circular convolution, with no hints as to the details for now $a[u]$ and $b[u]$ must be exactly evaluated.


Put another way - AFAIK, when we see equation-1, we simply assume that $a[u]$ and $b[u]$ are evaluated from $0 \leq u \leq N-1$, and then we do a circular convolution, modulo-$N$ - how/why then do we derive the proper evaluations for $a$ and $b$ from this equation??


Thanks!



Answer




You're interpretation that a[u] and b[u] are both N length sequences is incorrect. a[u] is of length N because it is effectively truncated by the length of your input sequence. The period of b[u] is not necessarily N. Depending on the value of N, b[u] may not even be periodic, so what you have in Equation 7 is a linear convolution of two sequences. One sequence (a[u] has length N, the other (b[u]) has length unknown. The limits of summation in the convolution equation are strictly set by the length of $a[u]$ and don't imply circular convolution or anything about the length of sequence b[u].


Next you have to consider what range of values of b[u] are required to compute the convolution. Straight from the form of EQUATION 1, you see that b[u] in general must be defined in the range -[N-1] to N-1. The negative boundary is set by the most negative b[u] value required to calculate X[0] (i.e. you need b[0−(N-1)] for this calculation). The positive boundary is set by the maximum spectral step you want to evaluate in the chirp Z transform (X[K] max). From the way the question has been posed, X[k] is best interpreted as $X[N-1]$.


The reason your computations work with a range that only goes back as far as -[M-1] is because you started with an M length sequence and zero padded it up to length N which zeros out any b[u] values going back further than -[M-1].


There is an interesting write up of Chirp-Z along with Goertzel's algorithm that I think puts this in proper perspective: http://www.google.com/url?sa=t&rct=j&q=goertzel%20algorithm&source=web&cd=3&sqi=2&ved=0CEIQFjAC&url=http%3A%2F%2Focw.mit.edu%2Fcourses%2Felectrical-engineering-and-computer-science%2F6-341-discrete-time-signal-processing-fall-2005%2Flecture-notes%2Flec20.pdf&ei=rawjUd6sM-LE0QHJ2IGQAQ&usg=AFQjCNEFEo7R7VSmk5jWYon0f4WqZYpvzg&bvm=bv.42553238,d.dmQ


No comments:

Post a Comment

readings - Appending 内 to a company name is read ない or うち?

For example, if I say マイクロソフト内のパートナーシップは強いです, is the 内 here read as うち or ない? Answer 「内」 in the form: 「Proper Noun + 内」 is always read 「ない...