## Convolution of Short Signals

Figure 7.1 illustrates the conceptual operation of filtering an input signal by a filter with impulse-response to produce an output signal . By the convolution theorem for DTFTs (§2.3.5),

*convolution*of and is defined by

In practice, we always use the DFT (preferably the FFT) in place of the DTFT, in which case we may write

*cyclic*(also called

*circular*) convolution [248]:

where means `` modulo .''

Another way to look at convolution is as the inner product of , and
, where
, *i.e.*,

*graphical convolution*in which the output sample at time is computed as an inner product of the impulse response after flipping it about time 0 and shifting time 0 to time . See [248, p. 105] for an illustration of graphical convolution.

### Cyclic FFT Convolution

Thanks to the convolution theorem, we have two alternate ways to perform cyclic convolution in practice:

- Direct calculation in the time domain using (7.2)
- Frequency-domain convolution:
- Fourier Transform both signals
- Perform term by term multiplication of the transformed signals
- Inverse transform the result to get back to the time domain

### Acyclic FFT Convolution

If we add enough trailing zeros to the signals being convolved, we can
obtain *acyclic* convolution embedded within a cyclic
convolution. How many zeros do we need to add? Suppose the signal
consists of contiguous nonzero samples at times **0** to
, preceded and followed by zeros, and suppose is nonzero
only over a block of samples starting at time 0. Then the
acyclic convolution of with reduces to

The number is easily checked for signals of length 1 since , where is 1 at time zero and 0 at all other times. Similarly,

When or is infinity, the convolution result can be as
small as 1. For example, consider
, with
, and
. Then
. This is an example of what is called *deconvolution*.
In the frequency domain, deconvolution always involves a pole-zero
cancellation. Therefore, it is only possible when or is
infinite. In practice, deconvolution can sometimes be accomplished
approximately, particularly within narrow frequency bands
[112].

We thus conclude that, to embed acyclic convolution within a cyclic
convolution (as provided by the FFT), we need to *zero-pad* both
operands out to length , where is at least the sum of the
operand lengths (minus one).

#### Acyclic Convolution in Matlab or Octave

In Matlab or Octave, the `conv` function implements acyclic
convolution:

octave:1> conv([1 2],[3 4]) ans = 3 10 8Note that it returns an output vector which is long enough to accommodate the entire result of the convolution, unlike the

`filter`primitive, which always returns an output signal equal in length to the input signal:

octave:2> filter([1 2],1,[3 4]) ans = 3 10 octave:3> filter([1 2],1,[3 4 0]) ans = 3 10 8

#### Pictorial View of Acyclic Convolution

Figure 7.2 shows schematically the result of convolving
two zero-padded signals and . In this case, the signal
starts some time after , say at . Since begins at
time **0**, the output starts promptly at time , but it takes some
time to ``ramp up'' to full amplitude. (This is the *transient
response* of the FIR filter .) If the length of is , then
the transient response is finished at time
. Next, when
the input signal goes to zero at time , the output reaches
zero samples later (after the filter ``decay time''), or time
. Thus, the total number of nonzero output samples is
.

If we don't add enough zeros, some of our convolution terms ``wrap
around'' and add back upon others (due to modulo indexing). This can
be called *time domain aliasing*. Zero-padding in the time
domain results in more samples (closer spacing) in the frequency
domain, *i.e.*, a higher `sampling rate' in the frequency domain. If we
have a high enough spectral sampling rate, we can avoid time aliasing.

The motivation for implementing acyclic convolution using a zero-padded cyclic convolution is that we can use the Fast Fourier Transform (FFT) to implement cyclic convolution when its length is a power of 2.

### Acyclic FFT Convolution in Matlab or Octave

The following example illustrates the implementation of acyclic convolution using the FFT in Matlab or Octave:

x = [1 2 3 4]; h = [1 1 1]; nx = length(x); nh = length(h); nfft = 2^nextpow2(nx+nh-1) xzp = [x, zeros(1,nfft-nx)]; hzp = [h, zeros(1,nfft-nh)]; X = fft(xzp); H = fft(hzp); Y = H .* X; format bank; y = real(ifft(Y)) % zero-padded result yt = y(1:nx+nh-1) % trim and print yc = conv(x,h) % for comparisonProgram output:

nfft = 8 y = 1.00 3.00 6.00 9.00 7.00 4.00 0.00 0.00 yt = 1.00 3.00 6.00 9.00 7.00 4.00 yc = 1 3 6 9 7 4

### FFT versus Direct Convolution

The following table compares the number of operations needed to perform the convolution of two length sequences for various values of :

N | FFT | Direct Convolution |

4 | 176 | 16 |

32 | 2560 | 1024 |

64 | 5888 | 4096 |

128 | 13,312 | 16,384 |

256 | 29,696 | 65,536 |

2048 | 311,296 | 4,194,304 |

In this example (adapted from [261]), the FFT (software) beats
direct time-domain convolution at length 128 and higher. It takes
approximately multiply/add operations to calculate the
convolution summation directly, while it takes on the order of
** log** operations to use the FFT method.
(Note, by the way, that
can be calculated once in
advance for time-invariant filtering operations.)

In digital audio, FIR filters are often hundreds of taps long. For such filters, the FFT method is much faster than direct convolution in the time domain.

#### Audio FIR Filters

FIR filters shorter than the ear's ``integration time'' can generally
be characterized by their magnitude frequency response (no perceivable
``delay effects''). The nominal ``integration time'' of the ear can
be defined as the reciprocal of a critical bandwidth of hearing.
Using Zwicker's definition of critical bandwidth
[282], the smallest critical bandwidth of hearing
is approximately 100 Hz (below 500 Hz). Thus, the nominal integration
time of the ear is 10ms below 500 Hz. (Using the
equivalent-rectangular-bandwidth (ERB) definition of critical
bandwidth [169,252], longer values are obtained).
At a 50 kHz sampling rate, this is 500 samples. Therefore, FIR
filters shorter than the ear's ``integration time,'' *i.e.*, perceptually
``instantaneous,'' can easily be hundreds of taps long (as discussed
in the next section). FFT convolution is consequently an important
implementation tool for FIR filters in digital audio applications.

#### Example 1: Low-Pass Filtering by FFT Convolution

In this example, we design and implement a length FIR low-pass filter having a cut-off frequency at Hz. The filter is tested on an input signal consisting of a sum of sinusoidal components at frequencies Hz. We'll filter a single input frame of length , which allows the FFT to be samples (no wasted zero-padding).

% Signal parameters: f = [ 440 880 1000 2000 ]; % frequencies M = 256; % signal length Fs = 5000; % sampling rate % Generate a signal by adding up sinusoids: x = zeros(1,M); % pre-allocate 'accumulator' n = 0:(M-1); % discrete-time grid for fk = f; x = x + sin(2*pi*n*fk/Fs); end

A plot of the synthesized input signal is shown in Fig.7.3. Next we design the lowpass filter using the window method:

% Filter parameters: L = 257; % filter length fc = 600; % cutoff frequency % Design the filter using the window method: hsupp = (-(L-1)/2:(L-1)/2); hideal = (2*fc/Fs)*sinc(2*fc*hsupp/Fs); h = hamming(L)' .* hideal; % h is our filter

Figure 7.4 plots the impulse response and amplitude response of our FIR filter designed by the window method. Next, the signal frame and filter impulse response are zero padded out to the FFT size and transformed:

% Choose the next power of 2 greater than L+M-1 Nfft = 2^(ceil(log2(L+M-1))); % or 2^nextpow2(L+M-1) % Zero pad the signal and impulse response: xzp = [ x zeros(1,Nfft-M) ]; hzp = [ h zeros(1,Nfft-L) ]; % Transform the signal and the filter: X = fft(xzp); H = fft(hzp);

Figure 7.5 shows the input signal spectrum and the filter amplitude response overlaid. We see that only one sinusoidal component falls within the passband.

Now we perform cyclic convolution in the time domain using pointwise multiplication in the frequency domain:

Y = X .* H;The modified spectrum is shown in Fig.7.6.

The final acyclic convolution is the inverse transform of the pointwise product in the frequency domain. The imaginary part is not quite zero as it should be due to finite numerical precision:

y = ifft(Y); relrmserr = norm(imag(y))/norm(y) % check... should be zero y = real(y);

Figure 7.7 shows the filter output signal in the time
domain. As expected, it looks like a pure tone in steady state. Note
the equal amounts of ``pre-ringing'' and ``post-ringing'' due to the
use of a linear-phase FIR filter.^{8.1}

For an input signal approximately samples long, this example is
2-3 times faster than the `conv` function in Matlab (which is
precompiled C code implementing time-domain convolution).

#### Example 2: Time Domain Aliasing

Figure 7.8 shows the effect of insufficient zero padding, which can
be thought of as *undersampling in the frequency domain*. We
will see *aliasing* in the time domain results.

The lowpass filter length is and the input signal consists of
an impulse at times and
, where the data frame
length is . To avoid time aliasing (*i.e.*, to implement
acyclic convolution using the FFT), we must use an FFT size at
least as large as
. In the figure, the FFT sizes ,
, and are used. Thus, the first case is heavily time
aliased, the second only slightly time aliased (involving only some of
the filter's ``ringing'' after the second pulse), and the third is
free of time aliasing altogether.

**Next Section:**

Convolving with Long Signals

**Previous Section:**

Summary