A Fast Guaranteed-Stable Sliding DFT Algorithm
This blog presents a most computationally-efficient guaranteed-stable real-time sliding discrete Fourier transform (SDFT) algorithm. The phrase “real-time” means the network computes one spectral output sample, equal to a single-bin output of an N‑point discrete Fourier transform (DFT), for each input signal sample.
Proposed Guaranteed Stable SDFT
My proposed guaranteed stable SDFT, whose development is given in [1], is shown in Figure 1(a). The output sequence Xk(n) is an N-point DFT's kth-bin spectral output sequence, where frequency index variable k is an integer in the range 0 ≤ k ≤ N–1, computed on a sample-by-sample basis. (The proposed Figure 1(a) SDFT is more computationally efficient than my guaranteed stable SDFT network presented in [1].)
Instead of using the traditional SDFT’s (see Appendix A in the PDF file associated with this blog) marginally stable complex resonator here I used a 2nd-order real resonator shown in Figure 1(a).
In Figure 1(b) we see the proposed network's real resonator has conjugate poles lying on the z‑plane's unit circle, for and N = 8 and k = 1, where the number 2 near the lower pole indicates dual zeros at $z = e^{-j2{\pi}k/N} = e^{-j{\pi}/4}$. The second of those dual zeros is produced by the feedforward stage of the network.
Due to pole/zero cancellation Figure 1(c) shows the network's equivalent z‑plane pole/zero locations. The proposed SDFT network's frequency magnitude response is given in Figure 1(d) for N = 8 and k = 1 where the peak response is located at $2{\pi}k/N = {\pi}/4$ radians/sample (corresponding to a cyclic frequency of $kf_s/N = f_s/8$ where fs is the x(n) input’s sample rate measured in Hz).
The Figure 1(a) network's time-domain impulse response is exactly k cycles of a complex exponential, whose frequency is $2{\pi}k/N$ radians/sample, having a duration of an N‑samples.
The Figure 1(a) network's transfer function is:$$H(z) = {e^{j2{\pi}k/N}-z^{-1}-e^{j2{\pi}k/N}z^{-N}+z^{-N-1}\over 1-2cos(2{\pi}k/N)z^{-1}+z^{-2}}.\qquad(1)$$The MATLAB/Octave code to display the frequency magnitude response of the Figure 1(a) network (for N = 8 and k = 1) is:
N = 8; k = 1;
B = [exp(j*2*pi*k/N), -1, zeros(1,N-2), -exp(j*2*pi*k/N), 1];
A = [1, -2*cos(2*pi*k/N), 1];
freqz(B, A, 512, 'whole')
Quantization of the center coefficient in the denominator in Eq. (1) will cause the real resonator's poles to shift, as will be discussed later, slightly from their desired positions along the unit circle. However, what is not widely known is that regardless of coefficient quantization the poles will always lie exactly on the z-plane’s unit circle guaranteeing resonator numerical stability [2]. So the proposed Figure 1 network produces the same spectral output samples as a single bin of an N‑point DFT with additional advantage that its output spectral samples are updated for each new input signal sample.
Proposed SDFT Computational Efficiency
The proposed Figure 1(a) guaranteed stable SDFT network is computationally more efficient than previously published guaranteed stable SDFT networks. That comparison is presented in Table 1.
The downloadable PDF file associated with this blog contains the following additional SDFT information:
• The effects of quantized coefficients on the Figure 1(a) SDFT,
• An SDFT signal tracking example,
• A brief tutorial on the traditional (unstable) SDFT,
• How to implement time-domain windowing with SDFTs,
• How to implement an SDFT for noninteger values of k.
Conclusions
I have proposed the most computationally efficient sliding DFT network, shown in Figure 1(a). That network computes a standard N‑point DFT's kth bin output on a real-time sample for sample basis and is guaranteed stable for all coefficient quantization scenarios.
References
[1] R. Lyons and C. Howard, “Improvements to the Sliding DFT Algorithm,” IEEE Signal Processing Mag., vol. 38, no. 4, pp. 119–127, July. 2021.
[2] R. Lyons, "Correcting an Important Goertzel Filter Misconception, dsprelated.com blog at: https://www.dsprelated.com/showarticle/796.php
[3] K. Duda, "Accurate, Guaranteed Stable, Sliding Discrete Fourier Transform,” IEEE Signal Processing Mag., vol. 27, no. 6, pp. 124–127, Nov. 2010.
[4] C. Orallo and I. Carugati, "Single Bin Sliding Discrete Fourier Transform", Chapter 2 of Fourier Transforms: High-tech Application and Current Trends, IntechOpen Limited, London, pp. 25-42, 2017.
- Comments
- Write a Comment Select to add a comment
Dear Rick,
thanks for your great blog article, it's very interesting! Especially the "doubling" trick.
In my conclusion, reducing the comb delay line length to N/2 in combination with doubling the frequency bins should result in the originally desired number of frequency bins.
This could be particularly advantageous in real-time audio applications, for example, where the shortest possible latency (impacted by the length of the comb delay line) is an important factor. But also frequency resolution as the counterpart to time resolution...
However, I'm not sure about the general "meaning" behind SDFT bin doubling.
Is it comparable to the classic FFT zero padding, where the increased number of resulting frequency bins is an effect of the interpolation and thus does not provide any higher frequency resolution than the original number of frequency bins?
Or would the doubling of the SDFT frequency bins lead to a higher frequency resolution, especially without the need for a larger comb delay line?
Juergen
Hello jurihock (Juergen).
Allow me to illustrate what I call "frequency doubling" with an example:
Let's consider a sliding DFT (SDFT) defined by N = 16 and the SDFT's input sequence is real-valued samples at a sample rate of 32 kHz. So we're only interested in the DFT resonant frequencies from zero Hz to Fs/2 Hz.
When the output of the comb filter's delay line is subtracted from the input sequence the available DFT bin-center frequencies are integer multiples of Fs/N (2000 Hz). Those frequencies are shown in the left column of the below table.
Now in an alternate implementation, when the output of the comb filter's
delay line is added to the input sequence the available DFT bin-center frequencies
are now the frequencies shown in
the right column of the above table.
So ..., for real-valued inputs, over the frequency range of zero to Fs/2 Hz we can set our SDFT bin center to be any of the frequencies in the left or right column of the above table. (No matter what SDFT bin-center frequency you choose, the SDFT's frequency magnitude response mainlobe's null-to-null width is 4000 Hz.)
Juergen, if you have any other questions for me please let me know.
Hi Rick,
When I use your Octave code for freqz:
N = 8; k = 1;
B = [exp(j*2*pi*k/N), -1, zeros(1,N-2), -exp(j*2*pi*k/N), 1];
A = [1, -2*cos(2*pi*k/N), 1];
freqz(B, A, 512, 'whole')
I get the following plot, am I doing it wrong? as it is different from your plot.
Hi kaz.
Your plot is correct. Sorry for the confusion. I plotted
the spectrum in my blog's Figure 1(d) with a linear
vertical axis and zero Hz in the middle of the x-axis.
I did that in hopes of improving the plot's simplicity
and clarity.
However, in my MATLAB code I used the standard standalone
'freqz()" command which plots the spectrum with a logarithmic
vertical axis and zero Hz is at the very left end of the
x-axis. Again, sorry for the confusion.
kaz, try out the following MATLAB code:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%N = 8; k = 1;
B = [exp(j*2*pi*k/N), -1, zeros(1,N-2), -exp(j*2*pi*k/N), 1];
A = [1, -2*cos(2*pi*k/N), 1];
[Spec, W] = freqz(B, A, 256, 'whole'); % Linear spectral vert. axis
Spec(k*256/N+1) = N; % Correct for the unwanted zero spec sample
Spec = fftshift(Spec); % Put zero Hz in the center of the x-axis
figure(1)
plot(linspace(-N/2, N/2, 256), abs(Spec))
ylabel('Linear'), xlabel('k'), grid on
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Regarding the above oddball 'Spec(k*256/N+1) = N;' command, try
inserting a '%' at the beginning of that command
(%Spec(k*256/N+1) = N;) and run the code to see why I used that command.
at a frequency exactly equal to 'k' it ends up trying to compute zero
divided by zero, which is "mathematically undefined". In which case
MATLAB sets that "undefined" spectral sample value equal to zero
So my 'Spec(k*256/N+1) = N;' command eliminates that plotting problem.
I give a detailed explanation of this "mathematically undefined"
spectral sample topic in Section 3.13.1 of my "Understanding DSP" textbook.
Thanks Rick,
so you are saying that freqz arrives at 0/0 and sets it to 0 but isn't 0/0 = 0?
So do you see that the function is failing? if so why it doesn't replace it like you did?
I know log of zero is undefined. But this is something else.
Hi kaz. Your reply is interesting.
I don't know what the ratio of 0/0 means but I'm pretty that 0/0 does not equal zero. The ratio of 0/0 seems to be one of those simple mathematical expressions that no one knows what it means. So they call it "undefined', or "indeterminate".
I have not tried to examine the internal code in MATLAB's 'freqz()' function. But I'm willing to wager a pint of Kentucky bourbon that when the internal code encounters the calculation sin(0)/sin(0) = 0/0 the code assigns that ratio to be equal to zero.
In my SDFT analysis I know to set 0/0 equal to 'N' but MATLAB's 'freqz()' function wouldn't know what to do with 0/0.
Hi kaz.
If I win the pint of bourbon mentioned in my previous reply, I'll bring the pint over to your place and we can share it. However, if you live in Australia then you'll have to pay for my airfare to travel to your place.
Thanks Rick,
That is kind of you, you can use Amazon delivery service to UK to avoid the hassle of travel and queues.
Back to that pesky zero and apologies if I am off track:
Matlab functions (some) seem to detect zero division and manage it for example, it sets result to 1 here if x is zero:
%%%%%%%%%%%%%%%%%%%%%%%%%%%
M = 39; Fc = 20; Fs = 100;
x = linspace(-(M-1)/2,(M-1)/2,M);
h = sinc(x*2*Fc/Fs);
%%%%%%%%%%%%%%%%%%%%%%%%%%
but if I go my way, I need to detect and set the value:
x = linspace(-(M-1)/2,(M-1)/2,M);
x(x==0) = 10^-30; %avoid division by zero
h = sin(x*2*pi*Fc/Fs)./(x*2*pi*Fc/Fs);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Detecting division by zero has to be followed by a replacement but requires knowledge as what to replace with. The best way is to replace at source with some small value as I did above and let the calculation ran over it. That is if we can.
However with 0/0 and it might sound funny: We all agree N/0 => inf, to me this is not true always? I am more happy with 0/0 = 1
proof: 10^-30/10^-30 = 1 and we can approach 0/0 in the limit
Imagine what what else would I say if I drink that pint !
by the way if you set freqz to odd resolution the problem is over.
Hi kaz.
I like your 'sinc()' and sin(x)/x examples.
MATLAB's built-in 'sinc(x)' function checks to see if 'x' is ever equal to zero, in which case it sets the function's output sample equal to 1.
So you're happy with 0/0 = 1, huh? If I had to pick a value for 0/0 I'd follow the guidance of the famous British novelist Douglas Adams who claimed the most important number in the universe is 42. So I say:
0/0 = 42 (1)
I verify the above Eq. (1) by multiplying both sides by zero giving us:
0 = 0*42
which is correct. :-)
Dear Rick,
Thanks for sharing, this is great! However, if I use this for bin 0 (DC) and send a single sample impulse into it, the ouput won't ever reach 0 again. Since cos(0) = 1 the resonator part will continue forever? After N samples it will receive the negative impulse: x(N) - x(N-N) = -1. From there only zeros follow but the resonator part won't decay, it remains constant.
Anything wrong with this pseudo code?
comb = x[n] - x[n - N]
resonator = comb + (2 * cos(0) * z1) - z2
z2 = z1
z1 = resonator
Of course z1 and z2 are initialised to 0.
Thank you!
Kind regards,
Gustav
Hell gustl.
If you implement my blog's Figure 1 network properly, if N = 20 and k = 0 then its impulse response will be 20 unity-valued samples followed by a sequence of all zero-valued samples. That impulse response is equivalent to a simple 20-point moving averager.
If you have an application that requires 'k' to be k = 0, then it would smarter to use a moving average filter rather than my Figure 1 SDFT network.
gustl, I'm sorry but I don't understand your pseudo code, so I am not able to comment on it.
Hi Rick,
Yes, that‘s right, nothing wrong about your proposed network. However, if you just look at the resonator part when using k=0 I‘m a bit worried that it could overflow at some point. But on the other hand it doesn‘t really matter because for k=0 a moving average will work fine, as you said.
By the way, I‘ve got an implementation of a sliding DFT using lookup tables. It needs only 2 multiplies and 4 adds per bin but it needs to be phase corrected which requires another multiply and add. If you are interested I can send it to you.
Kind regards,
Gustav
Hello Gustav (gustl).
Regarding your Reply's first paragraph, you are correct. My Figure 1 network has poor pole/zero cancellation at z = 1 when k = 0. But as we agree, for k = 0 a simple moving average is the smarter implementation.
I'd be happy to look at your sliding DFT using a lookup table. You're welcome send information regarding your sliding DFT using a lookup table to me at:
R.Lyonws@ieee.org.
Thanks.
Hi Rick,
Thanks for the interesting article. I am having trouble understanding two implementation details with your network above.
1. Would this network work for a complex-valued input signal? What would change if so? Would the real resonator portion only apply to the real part of the comb output, something like below?
2. In your reference [1] you discuss using non-integer values of k by adding a complex multiplication in the comb stage (Figure 5a). Would you be able to perform the same trick in this network? This also leads to the question above as the output of the comb will be complex-valued once this is implemented.
Thank you!
Adam Troyer
Hello Adam (amtroyer). I'm glad my blog is of some interest to you.
[Question# 1] You wrote, "Would this network work for a complex-valued input signal?" If your words "this network" mean my blog's Figure 1 network then the answer is yes. My blog's Figure 1 network does work with complex-valued inputs. The change from real-valued inputs to complex-valued inputs would be an increase in the required number of real-valued additions and multiplications per input sample.
I modeled your posted Reply's network and as far as I can tell it will not work properly with real-valued or complex-valued input signals.
[Question# 2] You wrote, "Would you be able to perform the same trick in this network?" Again, I don't know what your words "this network" mean? We're dealing with three networks here: (1) my blog's Figure 1 network, (2) your posted Reply's network, and (3) reference [1]'s Figure 5(a) network. So I don't know which of those networks your words "this network" is referring to.
Assuming your "this network" means my blog's Figure 1 network, the answer to your Question# 2 is, if you include a complex-multiply in my blog Figure 1 network's comb section then you can use a noninteger value for the 'k' variable.
Adam, if you have any further questions for me please let me know.
Hi Rick,
Sorry for the confusion with "this network". I should have chosen my words better! You were correct that I meant blog Figure 1 network in all of those instances.
Regarding question 1, I found the error I was making when using a complex-valued input and I have the blog Figure 1 network working correctly now for real and complex inputs.
Fixing this has brought me two more questions, if you don't mind my asking:
1. I am finding that with a complex-valued input the magnitude of the SDFT output is 2x larger than with a real-valued input. For example, if my input signal is a cos with amplitude of 1, then |X(k)|/(N/2) = 1 (after the initial N samples). If my input signal is a complex exponential with a magnitude of 1, |X(k)|/(N/2) = 2. Is this correct? If so, is this occurring because for a real-valued input half of the input signal power is in the -k bin, so it is not contributing to the output, while for a complex-valued input the entirety of the input signal power is in the k bin?
2. Why does the subtraction in the feedforward stage only get applied to the cos() leg and not the sin() leg?
Thanks for the help!
Adam
Hello Adam (amtroyer).
Regarding your first question: Your SDFT results are correct for a complex-valued input sequence. (That’s good news because it means your software modeling is correct. 😊) And your justification for your results is “essentially” correct. For a real-valued input sinusoid, containing exactly k cycles, half the spectral energy will be in the kth bin and the other half of the spectral energy will be in the (N-k)th bin.
As for your second question: I don’t have an answer, off the top of my head, on why there’s a subtraction in the feedforward stage’s real path and not in the imaginary path. The full details of how I derived my blog’s Figure 1 network are given in my blog’s Reference [1].
Hello Rick,
your article is really interesting (as every I had read to this moment - I am glad that you share your work in way that is understandable for most of us :) )
I have a question also connected to complex-valued SDFT solution based on Figure 1a network. (I have just read article and not yet checked this implementation and possibility to use it in algorithm which I ask below).
Is it able to use mentioned previously complex-valued SDFT in connection with algorithms that you proposed in article "Four ways to Compute an Inverse FFT Using the Forward FFT Algorithm" to implement SIDFT (Sliding Inverse Discrete Fourier Transform) or is there a better way to compute inverse transformation using proposed (mentioned in this article) SDFT algorithm?
Best regards,
Michal (Michael) Kaczmarek
Hello Michael. I just now saw your Reply here.
I've never thought about implementing a Sliding Inverse Discrete Fourier Transform (SIDFT) before now. But if I had to implement an SDFT the first thing I would try is the following:
Start with the following inverse DFT network.
Replace the above "Forward FFT" block with the following forward Sliding DFT network:
All of the processing in the forward SDFT network is complex-valued, so there are two comb filters, two 2nd-order resonators, six multipliers, etc.
Compute the real-valued "real(P(n))-imag(Q(n))" and "imag(P(n)) + real(Q(n))" sequences. Perform the final swapping shown in the inverse DFT network. And divide both real-valued swapped sequences by N.
Hello Rick,
Thank you very much for your response.
This was exactly my initial though on this topic. I will try this solution on my algorithm.
Best regards,
Michael
To post reply to a comment, click on the 'reply' button attached to each comment. To post a new comment (not a reply to a comment) check out the 'Write a Comment' tab at the top of the comments.
Please login (on the right) if you already have an account on this platform.
Otherwise, please use this form to register (free) an join one of the largest online community for Electrical/Embedded/DSP/FPGA/ML engineers:
