3
$\begingroup$

I have a discrete audio stream $x$ that needs to be processed in real-time. Specifically, as the each new sample is received, I would like to compute a Fourier transform of the last $n$ samples of the signal. So, for example, upon receiving the $i^{th}$ sample, I would like to compute the Fourier transform of the sequence $(x_{i-n+1}, x_{i-n+2}, \dots, x_{i-1}, x_{i})$. However this means that when I receive the $(i+1)^{th}$ sample, I will be computing the Fourier transform of $(x_{i-n+2}, x_{i-n+3}, \dots, x_{i}, x_{i+1})$, which is almost identical to the first sequence. So my question is, having already computed the Fourier transform of the first sequence, is there some way to compute the new Fourier transform from that, instead of having to recompute it from scratch each time?

$\endgroup$
  • 1
    $\begingroup$ regarding a practical answer, two questions arise, and I think they might both be very relevant when giving you an answer that actually helps you rather than just answer your question: 1. Why do you need each consecutive possible DFT? What is the application? 2. Please give us an idea of the sampling rates and DFT lengths. Orders of magnitude are sufficient. (the 1. question is way more important) $\endgroup$ – Marcus Müller Jan 26 '17 at 20:33
  • $\begingroup$ Normally, audio stream samples do not come into a processing system in "real-time", but in buffers or blocks of some length (for example 256 new samples every 6 milliseconds, or so, which is a non-zero latency). Thus, a full FFT will usually be faster than 256 steps of a sliding DFT algorithm. $\endgroup$ – hotpaw2 Jan 26 '17 at 20:49
  • $\begingroup$ Here is a simple algorithm I will provide which would improve the computational efficiency by $\log_2(N)$ times where N is the length of the window whose contents are shifted to left by 1 after arrival of each new new sample to the rightmost position. for N = 1024 this would mean a boost of almost 10x compared to a direct N-point FFT which would require $N \log_2(N)$ complex MACs in the rough. $\endgroup$ – Fat32 Jan 26 '17 at 21:52
  • 1
    $\begingroup$ @MarcusMüller You're right, this isn't exactly the problem that I'm trying to solve. Basically, my goal is to write a real-time VAD algorithm which is based on computing windowed cepstra of the incoming audio stream. So I'm asking about online variants of the Fourier transform because that's one of the steps that's needed to compute the cepstrum. However it's true that I would rather not have to recompute this cepstrum at every new sample. Ideally, I would be computing the cepstrum over 40ms windows with 10ms overlaps with the signal being sampled at a rate of 16kHz. $\endgroup$ – jon_simon Jan 26 '17 at 22:44
  • 2
    $\begingroup$ Maybe you'd want to accept one of the answers below, because they answer the question you've asked. Then, ask a new question about solving your actual problem; at a one quarter overlap of a 640-point FFT, it hardly pays off to do any shiftin/out trickery. $\endgroup$ – Marcus Müller Jan 26 '17 at 22:55
7
$\begingroup$

(Note: the paper pointed by hotpaw2's link is actually describing in more detail the algorthm I presented here)

Consider a data window length of $N$ samples from $n=0$ to $n=N-1$. Let your original data window be $x_1[n]$, whose first sample is $x_{old} = x_1[0]$. Now your new data set is denoted as $x_2[n]$ whose samples are actually one sample left shifted verisons of $x_1[n]$, i.e. $x_2[n] = x_1[n+1]$ for $n=0,1,...,N-2$ plus a new arrived data to position $n=N-1$ which is denoted as $x_{new} = x_2[N-1]$.

Then the following algorithm will compute the N-point DFT, $X_2[k]$ of the new data set $x_2[n]$ from that of the already computed and stored N-point DFT $X_1[k]$ of the old data set $x_1[n]$ as:

$$ X_2[k] = e^{j \frac{2\pi}{N}k }( X_1[k] + (x_{new} - x_{old}) ) $$

for each $k=0,1,...,N-1$. This updated computation of $X_2[k]$ from pre-computed $X_1[k]$ requires N complex multiplications and N real additions. Compared to a direct N-point FFT which requires $N \log_2(N)$ complex MACs, this would therefore be an improvement by a factor of $\log_2(N)$ which for example would render to a factor of 10 when $N=1024$ using low level languages such as C.

$\endgroup$
4
$\begingroup$

The algorithm for which you may be looking is called a "sliding DFT". For a small number of result bins, that number of "sliding Goertzel" filters can also be used. Here's one online description: https://www.dsprelated.com/showarticle/776.php of how to implement a sliding DFT.

Basically, for each new sample, you add a new twiddle-multiplied vector to the summation vector and subtract a historical (saved and/or recalculated) vector to remove the oldest slice of the summation, thus keeping the number of vector slices in the DFT summation the same. Numerical noise (quantization and rounding, etc.) can cause this process to be unstable over a sufficiently large number of steps, so an occasional fresh full DFT may be required to check for drift.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.