1
$\begingroup$

Assume I want to filter a stream of data with a rather small FIR filter (in the order of 16 to 64 coefficients). In any case I have to split the input stream into chunks (blocks) of data.

Searching the web immediately leads to the "fast convolution",i.e.: zero-padding in time domain-> doing the FFTs -> multiply -> doing the IFFT -> overlap add or overlap save to get the linear convolution.

My question for understanding is: Suppose I do not care about the better efficiency that I get by the above process using the frequency domain or I simply want to avoid doing an DFT: is it in principle possible to just calculate the (linear) convolution sum in the time domain between the blocks of input data and the FIR impulse response and then combining the results of the linear convolutions in a similar manner than in overlap-add or overlap-save? If yes, can someone provide a pseudocode/algorithm?

$\endgroup$
4
  • $\begingroup$ Overlap-add, overlap-save and time-domain FIR filtering all produce the same output. So, you are free to choose the one that suits your needs best. Time-domain FIR filtering is nothing more than linear convolution, which can be implemented a number of ways (most commonly "direct-form", followed by "transposed direct-form"). See, for example, here. $\endgroup$
    – Harry
    Commented Oct 22, 2023 at 15:17
  • $\begingroup$ @Harry: Maybe I was not specific enough: In case, I do the convolution in time domain, the output of each "block convolution" is longer than the input block. According to which algorithms do I combine the outputs of the single block convolutions to get the time domain output stream? Is it the very same than with overlap-add/save? $\endgroup$
    – Junius
    Commented Oct 22, 2023 at 16:26
  • $\begingroup$ @Junius: Let's say your block size is M and your filter length is N. Then each block you calculate exactly M output samples using M+N-1 input samples. That means you need to keep the last N-1 samples from the previous block around. That's called "state keeping". $\endgroup$
    – Hilmar
    Commented Oct 22, 2023 at 16:37
  • $\begingroup$ @Junius I think Hilmar read between the lines and saw that you're asking specifically about how to implement this in software. This wasn't clear to me. But, in that case, Hilmar is correct that you just need to manage (save and restore) the state of the filter's internal "delay line" between each batch of data. This is not the same as overlap add/save. $\endgroup$
    – Harry
    Commented Oct 22, 2023 at 18:05

1 Answer 1

1
$\begingroup$

is it in principle possible to just calculate the (linear) convolution sum in the time domain

Yes, of course. Any standard FIR implementation can do this. Depending on what language and/or platform you are using, you can easily find code on the Internet for that.

The actual core algorithm is really quite simple: it's a multiply/accumulate operation. If you flip the impulse response this becomes a simple dot product which many platforms, processors, and/or compilers have hardware acceleration and/or optimized libraries for

The more tricky part is the state keeping. If the length of the impulse response is $N$, then you need to keep the recent most $N-1$ input samples around from block to block. That state is generally not consecutive with your input data so there is some non-trivial memory management required. Options are

  1. Implement a circular buffer for the input
  2. Pre-pend the state in front of the input
  3. Use a double sized state buffer and split the convolution into two dot-products.

Again, the best choice depends a lot on your specific details and HW target. For example circular buffers can be extremely efficient or horribly inefficient depending on how they are implemented and what the hardware supports. Personally, I like option 3 as it's a good compromise between efficiency and portability.

Below is a piece of Matlab code that demonstrates this. Unfortunately Matlab starts index array at 1 (as compared to 9) so some of the indices have a +1 or -1 which isn't entirely intuitive.

%% FIR filter with double size state buffer
nFIR = 16; % fir filter size
nBlock = 256; % input block size
nx = 2^16;  % lentgh of input vector;

% create filter and impulse response

h = randn(nFIR,1);  % random coefficients
x = randn(nx,1);    % random signal
y = zeros(size(x)); % output buffer
numBlocks = nx/nBlock;

% setup the state
state = zeros(2*nFIR-1,1);  % twice th elnegth of the filter -1 
inputPointer = 0;
tBlock = (1:nBlock)';  % block sized index vector, starting at 1
tFIR = (1:nFIR)';      % FIR sized index vector, starting at 1
hFlipped = flip(h)';  % flip and transpose
% loop over all blocks
for ib = 1:nBlock
  x0 = x(inputPointer+tBlock,1); % get next block of input data
  y0 = zeros(size(x0));  % initialze output block
  % copy the first N input samples behind the state
  state(nFIR-1+tFIR) = x0(tFIR,1);
  % inner loop 1: over the double state buffer
  for i = 1:nFIR
    y0(i) = hFlipped*state(i+tFIR-1);  % this is the dot product
  end
  % inner loop 2: over the input data buffer
  for i = 1:nBlock-nFIR
    y0(i+nFIR) = hFlipped*x0(i+tFIR);  % dot product
  end
  % save output block
  y(inputPointer+tBlock,1) = y0;
  inputPointer = inputPointer + nBlock;
  % update the state, save the last nFIR-1 samples
  state(1:nFIR-1) = x0(nBlock-nFIR+2:nBlock,1);
  
end

%% compare to the reference
yRef = fftfilt(h,x);
err = y-yRef;
fprintf('Error = %6.2fdB\n',10*log10(sum(err.^2)./sum(yRef.^2)));
 
$\endgroup$
4
  • $\begingroup$ The provided code is more informative than using built-in functions. However, it's worth noting that MATLAB's built-in filter() function allows you to set an initial state zi and save a final state zf using this syntax: [y,zf] = filter(b,a,x,zi). It's for doing exactly this type of batch processing. So, if you actually want to implement this in MATLAB, then using filter() will give you a tidier and faster implementation. $\endgroup$
    – Harry
    Commented Oct 22, 2023 at 18:12
  • $\begingroup$ @Hilmar: This perfectly answers the question, although I still have to find the time to think about some of the code details to make sure I understand :) $\endgroup$
    – Junius
    Commented Oct 26, 2023 at 22:04
  • $\begingroup$ @Hilmar: when you say "That state is generally not consecutive with your input data" -> do you mean that the input block size is in general not the same length as the filter length or do you mean somthing different? $\endgroup$
    – Junius
    Commented Oct 27, 2023 at 13:33
  • $\begingroup$ @Junius: good question! What I mean is that state and input data are typically not consecutive in the internal memory of the processor. If they were, the implementation would be a LOT easier and more efficient. There are ways to achieve this but this depends a lot on the language and hardware platform. $\endgroup$
    – Hilmar
    Commented Oct 28, 2023 at 1:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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