Sunday, November 9, 2014

Linear phase FIR filters: 4 types

This post is about the 4 types of linear phase finite impulse response (FIR) filters. The possibility to have a causal filter with an exactly linear phase response can be an advantage of FIR filters over infinite impulse response (IIR) filters for certain applications.
  1. Type I: symmetrical impulse response, odd length
  2. Type II: symmetrical impulse response, even length
  3. Type III: asymmetrical impulse response, odd length
  4. Type IV: asymmetrical impulse response, even length
When choosing one of these 4 types of linear phase filters there are mainly 3 things to consider:
  1. constraints on the zeros of the transfer function $H(z)$ at $z=1$ (DC) and $z=-1$ (Nyquist)
  2.  integer/non-integer group delay
  3.  phase shift (apart from the linear phase)
For type I filters (odd number of taps, even symmetry) there are no constraints on the zeros at $z=1$ and $z=-1$, the phase shift is zero (apart from the linear phase), and the group delay is an integer value.

Type II filters (even number of taps, even symmetry) always have a zero at $z=-1$ (i.e., half the sampling frequency), they have a zero phase shift, and they have a non-integer group delay.

Type III filters (odd number of taps, odd symmetry) always have zeros at $z=1$ and $z=-1$ (i.e. at $f=0$ and $f=f_s/2$), they have a 90 degrees phase shift, and an integer group delay.

Type IV filters (even number of taps, odd symmetry) always have a zero at $z=1$, a phase shift of 90 degrees, and a non-integer group delay.

This implies (among other things) the following:
  • Type I filters are pretty universal, but they cannot be used whenever a 90 degrees phase shift is necessary, e.g. for differentiators or Hilbert transformers.
  • Type II filters would normally not be used for high pass or band stop filters, due to the zero at $z=-1$, i.e. at $f=f_s/2$. Neither can they be used for applications where a 90 degrees phase shift is necessary.
  • Type III filters cannot be used for standard frequency selective filters because in these cases the 90 degrees phase shift is usually undesirable. For Hilbert transformers, type III filters have a relatively bad magnitude approximation at very low and very high frequencies due to the zeros at $z=1$ and $z=-1$. On the other hand, a type III Hilbert transformer can be implemented more efficiently than a type IV Hilbert transformer because in this case every other tap is zero.
  • Type IV filters cannot be used for standard frequency selective filters, for the same reasons as type III filters. They are well suited for differentiators and Hilbert transformers, and their magnitude approximation is usually better because, unlike type III filters, they have no zero at $z=-1$.
  • In some applications an integer group delay is desirable. In these cases type I or type III filters are preferred.

Saturday, November 8, 2014

FFT: the Cooley-Tukey algorithm for general composite lengths

The Cooley-Tukey fast Fourier transform (FFT) algorithm is probably the most well-known FFT algorithm. Their paper "An algorithm for the machine calculation of complex Fourier series" was published in 1965, and it is said to have caused the boom of digital signal processing (DSP) in the 1970s. Note, however, that the algorithm was already known to Gauss at the beginning of the 19th century. The algorithm is based on recursively splitting the problem into smaller sub-problems, and thereby decreasing the computational complexity. In its most well-known form (the radix-2 FFT), the algorithm reduces the complexity from $N^2$ to the famous $N\log(N)$ for a discrete Fourier transform (DFT) of length $N=2^M$.

In this article I want to discuss the Cooley-Tukey algorithm for general composite $N=N_1N_2$. The easiest way to see how the algorithm works is to use index mapping. Consider the length-$N$ DFT of a sequence $x[n]$: $$X[k]=\sum_{n=0}^{N-1}x[n]W_N^{nk}\tag{1}$$ with $W_N=e^{-j2\pi/N}$. Let's assume that $N$ can be factored as $N=N_1N_2$. Now use the following index mapping $$n=n_1N_2+n_2,\quad n_1=0,1,\ldots, N_1-1,\quad n_2=0,1,\ldots,N_2-1\\ k=k_1+N_1k_2,\quad k_1=0,1,\ldots, N_1-1,\quad k_2=0,1,\ldots,N_2-1\tag{2}$$ Plugging (2) into (1) gives $$X[k_1+N_1k_2]=\sum_{n_2=0}^{N_2-1}\sum_{n_1=0}^{N_1-1}x[n_1N_2+n_2]W_N^{(n_1N_2+n_2)(k_1+k_2N_1)}\tag{3}$$ The complex exponentials in (3) can be simplified as follows: $$W_N^{(n_1N_2+n_2)(k_1+k_2N_1)}=W_N^{N_2n_1k_1}W_N^{N_1N_2n_1k_2}W_N^{n_2k_1}W_N^{N_1n_2k_2}=\\= W_{N_1}^{n_1k_1}\cdot 1\cdot W_{N}^{n_2k_1}W_{N_2}^{n_2k_2}\tag{4}$$ because $W_N^{N_1}=W_{N_2}$, $W_N^{N_2}=W_{N_1}$, and $W_N^{N_1N_2}=W_N^N=1$. Plugging (4) back into (3) and rearranging terms gives $$X[k_1+N_1k_2]=\sum_{n_2=0}^{N_2-1}\left[W_N^{n_2k_1}\left(\sum_{n_1=0}^{N_1-1}x[n_1N_2+n_2]W_{N_1}^{n_1k_1}\right)\right]W_{N_2}^{n_2k_2}\tag{5}$$ Equation (5) can be interpreted as follows: the result can be computed by computing $N_2$ length-$N_1$ DFTs of subsampled versions of the input sequence $x[n]$ (this is the inner sum in (5)), then multiplying the result by twiddle factors $W_N^{n_2k_1}$, and finally computing $N_1$ length-$N_2$ DFTs (the outer sum in (5)). This procedure is most easily visualized using matrices:
  1. read the data row-wise into a $N_1\times N_2$ matrix
  2. compute DFTs over all $N_2$ columns (this can be done in-place)
  3. multiply the result element-wise with twiddle factors
  4. compute DFTs over all $N_1$ rows (in-place)
  5. read out the data column-wise to get the final length-$N$ result
Let's see if and how this procedure reduces the complexity of computing a length-$N$ DFT. Direct computation of a length-$N$ DFT requires $N^2$ complex multiplications. With the above procedure we need $N_2$ DFTs of length $N_1$ ($N_2N_1^2$ complex multiplications), $N_1N_2$ complex multiplications for the twiddle factors, and $N_1$ DFTs of length $N_2$ ($N_1N_2^2$ complex multiplications), resulting in a total of $$N_2N_1^2+N_1N_2+N_1N_2^2 = N_1N_2(N_1+N_2+1) = \\ = N(N_1+N_2+1) < N^2,\quad\text{for } N_1,N_2 > 2$$ So the procedure results in reduced complexity compared to the direct computation of a length-$N$ DFT. If the factors $N_1$ and $N_2$ themselves can be factored, then the method can be used recursively, which will result in a greater reduction of complexity. E.g., if $N=2^M$ this results in the famous radix-2 version of the Cooley-Tukey algorithm, which is discussed in most DSP textbooks.

In practice, an FFT routine which implements $N=2^M$ will be sufficient most of the time, because in most cases one is free to choose the DFT length, and the data can usually be zero-padded to the next power of two. However, there are cases where a non-power of 2 length is required, such as in 3GPP Long Term Evolution (LTE), where for a transmission bandwidth of 15MHz an FFT length of 1536 is required. Suppose that on a certain platform you have an FFT routine for $N=2^M$. In this case, since $1536=3\times 512$, the representation (5) can be used to split up the computation of the DFT in $3$ FFTs of length $512$, which can be computed by the available FFT routine, and $512$ length-3 FFTs, the implementation of which is trivial.

This little octave/Matlab script illustrates the procedure for the case $N=1536$:

 N = 1536;         % non-power of 2 FFT length used in LTE  
 N1 = 512; N2 = 3; % 512 x 3 = 1536  
 x = randn(N,1);   % some input signal  
   
 % generate N1 x N2 matrix of input data, reading them row wise  
 X = reshape(x,[N2,N1]).';  
   
 % DFT over each column  
 X = fft(X);  
   
 % Multiplication with twiddle factors  
 X = X.*exp(-1i*2*pi/N*(0:N1-1)'*(0:N2-1));  
   
 % DFT over each row  
 X = (fft(X.')).';  
   
 % read out column wise  
 X = X(:);

So even if you don't want to program your own FFT algorithm it's very useful to understand the Cooley-Tukey method because it can be used to implement non-power-of-two FFTs using basic (power-of-2) FFT routines, if the FFT length $N$ is composed of only a few (small) factors (other than 2), and this will be the case for almost all practical applications.

Friday, November 7, 2014

Discrete-time signals: time shift vs phase shift

This article is based on my answer to this question on dsp.stackexchange.com. There was also another question asking basically the same thing, so it seems that people get confused about the relation between time shifts and phase shifts for discrete-time signals.

First let's consider a continuous time sinusoidal signal: $$x(t)=\cos(\omega_0t)$$ Applying a time shift of $t_0$ gives $$x(t+t_0)=\cos(\omega_0(t+t_0))=\cos(\omega_0t+\omega_0t_0)=\cos(\omega_0t+\phi)$$ with $\phi=\omega_0t_0$. So for a continuous-time sinusoidal signal a time shift always corresponds to a phase shift, but, more importantly, the opposite is also always true: a phase shift always corresponds to a time-shift. The time shift is given by
$$t_0=\frac{\phi}{\omega_0}$$ With discrete-time signals, things are different. A time-shift always corresponds to a phase shift, but the opposite is generally not true. Let $$x[n]=\cos(n\theta_0)$$ Applying a time-shift of $n_0\in\mathbb{Z}$ gives $$x[n+n_0]=\cos((n+n_0)\theta_0)=\cos(n\theta_0+n_0\theta_0)$$ which corresponds to a phase shift of $\phi=n_0\theta_0$. However, assume another signal $$y[n]=\cos(n\theta_0+\phi),\quad\phi\in\mathbb{R}$$ If $n$ were a continuous variable, the signal $y[n]$ could be obtained from $x[n]$ by applying an appropriate time-shift. However, in discrete-time $y[n]$ cannot generally be obtained by time-shifting  $x[n]$. This is only possible if $$\phi=n_0\theta_0$$ is satisfied for some integer $n_0$. This can be seen as follows: $$y[n]=\cos(n\theta_0+\phi)=\cos((n+\phi/\theta_0)\theta_0)$$ which only equals $x[n+n_0]$ if $n_0=\phi/\theta_0$ is an integer.

In sum, in continuous time a phase shift corresponds to a time shift and vice versa. In discrete-time, a time shift always corresponds to a phase shift, but the opposite is generally not true.

Digital Filter Design: An Overview

I wrote this little article as an answer to a question on dsp.stackexchange.com and I think it gives a good overview for someone who is just starting with the topic. Enjoy reading and leave comments if you wish.

Digital filter design is a very large and mature topic and there is a lot of material available. What I want to try here is to get a beginner started and to make the existing material more accessible. Instead of digital filters I should actually be talking about discrete-time filters because I will not consider coefficient and signal quantization here.

A non-recursive linear time-invariant (LTI) filter can be described by the following difference equation $$y(n)=h_0x(n)+h_1x(n-1)+\ldots +h_{N-1}x(n-N+1)=\\=\sum_{k=0}^{N-1}h_kx(n-k)\tag{1}$$ where $y(n)$ is the output sequence, $x(n)$ is the input sequence, $n$ is the time index, $h_k$ are the filter coefficients, and $N$ is the filter length (the number of taps). The filter taps $h_k$ are also the impulse response of the filter because if the input signal is an impulse, i.e. $x(n)=\delta(n)$, then $y(n)=h_n$ (if the filter's memory has been initialized with zeros). Equation (1) describes a linear time-invariant finite impulse response (FIR) system. The sum on the right-hand side of (1) is a convolution sum, i.e. the output signal is obtained by convolving the input signal with the impulse response. This is always true, but for IIR filters we cannot explicitly compute the convolution sum because the impulse response is infinitely long, i.e. there are infinitely many coefficients $h_k$. One important advantage of FIR filters is that they are always stable, i.e. for a bounded input sequence, the output sequence is always bounded. Another advantage is that FIR filters can always be realized with an exactly linear phase, i.e. they will not add any phase distortion apart from a pure delay. Furthermore, the design problem is usually easier, as we will see later.

A recursive LTI filter is described by the following difference equation: $$y(n)=b_0x(n)+b_1x(n-1)+\ldots+b_Mx(n-M)-\\ -a_1y(n-1)-\ldots-a_Ny(n-N)\tag{2}$$
Equation (2) shows that the output is not only composed of weighted and delayed input samples, but also of weighted past output samples. In general, the impulse response of such a system is infinitely long, i.e. the corresponding system is an IIR (infinite impulse response) system. However, there are special cases of recursive filters with a finite impulse response. Note that the impulse response is not anymore given by either the coefficients $b_k$ or $a_k$ as in the case of FIR filters. One advantage of IIR filters is that steep filters with high stopband attenuation can be realized with much fewer coefficients (and delays) than in the FIR case, i.e. they are computationally more efficient. However, one needs to be careful with the choice of the coefficients $a_k$ because IIR filter can be unstable, i.e. their output sequence can be unbounded, even with a bounded input sequence.

Filters can be designed according to specifications either in the time (sample) domain or in the frequency domain, or both. I'll restrict myself to specifications in the frequency domain. In this case one needs to consider the frequency responses of FIR and IIR systems. The frequency response of a system is the Fourier transform of its impulse response, assuming that it exists (which is the case for stable systems). The frequency response of an FIR filter is $$H(e^{j\theta})=\sum_{k=0}^{N-1}h_ke^{-jk\theta}\tag{3}$$ where $\theta$ is the discrete-time frequency variable: $$\theta=\frac{2\pi f}{f_s}$$ with the actual frequency $f$ and the sampling frequency $f_s$. From (3) you can see that approximating a desired frequency response by an FIR system is basically a problem of polynomial approximation. For recursive systems we have $$H(e^{j\theta})=\frac{\sum_{k=0}^Mb_ke^{-j\theta}}{1+\sum_{k=1}^Na_ke^{-j\theta}}\tag{4}$$ and you get a rational approximation problem, which is usually much more difficult than the polynomial approximation problem in the case of FIR filters. From (3) and (4) you can see that the frequency response of an FIR filter is of course only a special case of the response of a recursive filter with coefficients $a_k=0$, $k=1,\dots,N$.

Let's now take a quick look at filter design methods. For FIR filters you could take an inverse Fourier transform of the desired frequency response to get the impulse response of the filter, which directly corresponds to the filter coefficients. Since you approximate the desired response by a finite length impulse response you should apply a smooth window to the obtained impulse response to minimize oscillations in the actual frequency response due to Gibbs' phenomenon. This method is called frequency-sampling method. For simple standard filters like ideal lowpass, highpass, bandpass or bandstop filters (and a few others), you could even analytically calculate the exact impulse response by taking the inverse Fourier transform of the ideal desired response: $$h_k=\frac{1}{2\pi}\int_{-\pi}^{\pi}H(e^{j\theta})e^{jk\theta}d\theta$$ This integral is easy to evaluate for piecewise constant desired responses, as is the case for ideal frequency-selective filters. This will give you an infinitely long, non-causal impulse response, which needs to be windowed and shifted to make it finite and causal. This method is know as window-design. There are of course many other FIR filter design methods. One important numerical method is the famous Parks-McClellan exchange algorithm which designs optimal filters with constant passband and stopband ripples. It is a numerical approximation method and there are many software implementations available, e.g. in Matlab and Octave. The most common IIR design method for frequency selective filters is the bilinear transformation method. This method simply uses analytical formulas for the design of optimal analog filters (such as Butterworth, Chebyshev, Cauer/elliptic, and Bessel filters), and transforms them to the discrete-time domain by applying a bilinear transformation to the complex variable $s$ (analog domain) which maps the (imaginary) frequency axis of the complex $s$-plane to the unit circle in the complex $z$-plane (discrete-time domain). Don't worry if you do not yet know much about complex transfer functions in the analog or discrete-time domain because there are good implementations available of the bilinear transform method, e.g. in Matlab or Octave.

There are of course many more interesting and useful methods, depending on the type of specifications, but I hope that this will get you started and will make any material you come across more understandable. A very good (and free) book covering some basic filter design methods (and a lot more) is Introduction to Signal Processing by Orfanidis. You can find several design examples there. Another great classic book is Digital Filter Design by Parks and Burrus.