So I want to continuously estimate the frequency content of a signal using the FFT such that at every sample timestep, I add a new sample to my current set of samples and remove the oldest sample in that set and perform an FFT on the resulting new set of samples. I was wondering is there any way to make things more efficient by reusing calculations I did for the old set of samples for the very similar push-pop new sample set?
$\begingroup$
$\endgroup$
2
-
2$\begingroup$ There is something called the "sliding DFT". Similar to the moving average, there are numerical issues of subtracting exactly what was added when an old sample falls offa the edge and is replaced by the new sample coming in. Sometimes there is error that does not go away. A solution to that is to provide some kinda decay. This is all closely related to the notion of Truncated IIR filters. $\endgroup$– robert bristow-johnsonCommented Apr 2 at 15:42
-
$\begingroup$ Usually you want to apply a non-trivial window prior to taking the FFT. I imagine that this makes clever «FFT reuse schemes» difficult. $\endgroup$– Knut IngeCommented Apr 3 at 2:47
Add a comment
|
1 Answer
$\begingroup$
$\endgroup$
I was wondering is there any way to make things more efficient by reusing calculations I did for the old set of samples
As long as you just want to do straight FFT with a shift one sample per frame that's quite easy to do.
Let's say we have signal $x[n]$ of length $N$ and it's DFT $X[k]$. When we get a new sample $a$, we prepend this to the existing signal and drop the last one, i.e. we form a new signal
$$y[n] = \begin{cases} a & n = 0 \\ x[n-1] & 1 \leq n < N \end{cases} $$
Then we can calculate it's DFT
$$Y[k] = X[k]e^{-j2\pi\frac{k}{N}} + a - x[N-1] \tag{1}$$
Here is roughly why this works: Multiplication with $e^{-j2\pi\frac{k}{N}}$ is simply a circular shift of the time domain buffer by one sample. This rotates the last sample into the first position and the operation $a - x[N-1] $ replaces this with the new sample. Now we would need to take the FFT of this correction term but the FFT of a single values at $n=0$ is just a constant of the same value in the frequency domain.
Here is some Matlab code
%% FFT with one ssmple increment
nx = 1024; % number of samples
fDelay = exp(-2i*pi*(0:nx-1)'/nx); % one sample delay transfer functions
x0 = randn(nx,1); % random signal
fx0 = fft(x0); % and it's DFT
% create a new sample
a = randn(1,1);
% cacluated updsated FFT
fy = fx0.*fDelay + a - x0(nx);
%% Error check
y = [a; x0(1:nx-1)]; % prepend new sample and drop the last one
fref = fft(y); % FFT
err = fy - fref; % Difference
fprintf('Error = %6.2fdB\n',20*log10(rms(err)./rms(fref)));
