FIR Design

This section describes the finite impulse response filter design capabilities in liquid. This includes basic low-pass filter design using the windowed-\(\sinc\) method, square-root Nyquist filters, arbitrary design using the Parks-McClellan algorithm, and some useful miscellaneous functions.

Window prototype∞

The ideal low-pass filter has a rectangular response in the frequency domain and an infinite \(\sin(t)/t\) response in the time domain. Because all time-dependent filters must be causal, this type of filter is unrealizable; furthermore, truncating its response results in poor pass-band ripple stop-band rejection. An improvement over truncation is offered by use of a band-limiting window. Let the finite impulse response of a filter be defined as

$$ h(n) = h_i(n) w(n) $$

where \(w(n)\) is a time-limited symmetric window and\(h_i(n)\) is the impulse response of the ideal filter with a cutoff frequency \(\omega_c\) , viz.

$$ h_i(n) = \frac{\omega_c}{\pi} \left( \frac{\sin\omega_c n}{\omega_c n} \right), \,\,\,\, \forall n $$

A number of possible windows could be used; the Kaiser window is particularly common due to its systematic ability to trade transition bandwidth for stop-band rejection. The Kaiser window is defined as

$$ w(n) = \frac{ I_0\left[\pi\alpha\sqrt{1-\left(\frac{n}{N/2}\right)^2}\right] }{ I_0\left(\pi\alpha\right) } \,\,\,\, -N/2 \leq n \leq N/2, \,\, \alpha \geq 0 $$

where \(I_\nu(z)\) is the modified Bessel function of the first kind of order \(\nu\) and \(\alpha\) is a shape parameter controlling the window decay.\(I_\nu(z)\) can be expanded as

$$ I_\nu(z) = \left( \frac{z}{2} \right)^\nu \sum_{k=0}^{\infty}{ \frac{ \left(\frac{1}{4}z^2\right)^k }{ k!\Gamma(k+\nu+1) } } $$

The sum in [ref:eqn-filter-firdes-besseli_infinite_sum] converges quickly due to the denominator increasing rapidly, (and in particular for \(\nu=0\) the denominator reduces to \((k!)^2\) ) and thus only a few terms are necessary for sufficient approximation. The sum [ref:eqn-filter-firdes-besseli_infinite_sum] converges quickly due to the denominator increasing rapidly, thus only a few terms are necessary for sufficient approximation. For more approximations to \(I_0(z)\) and \(I_\nu(z)\) , see [ref:section-math] in the math module. Kaiser gives an approximation for the value of \(\alpha\) to give a particular sidelobe level for the window as {cite:Vaidyanathan:1993((3.2.7))}

$$ \alpha = \begin{cases} 0.1102 (A_s - 8.7) & A_s \gt 50 \\ 0.5842 (A_s - 21)^{0.4} & 21 \lt A_s \le 50 \\ 0 & \text{else} \end{cases} $$

where \(A_s \gt 0\) is the stop-band attenuation in decibels. This approximation is provided in liquid by the kaiser_beta_As() method, and the length of the filter can be approximated with estimate_req_filter_len() (see [ref:section-filter-misc] for more detail on these methods).

The entire design process is provided in liquid with the firdes_kaiser_window() method which can be invoked as follows:

liquid_firdes_kaiser(_n, _fc, _As, _mu, *_h);

where _n is the length of the filter (number of samples), _fc is the normalized cutoff frequency (\(0 \leq f_c \leq 0.5\) ), _As is the stop-band attenuation in dB (\(A_s \gt 0\) ), _mu is the fractional sample offset (\(-0.5 \leq \mu \leq 0.5\) ), and *_h is the \(n\) -sample output coefficient array. Listed below is an example of the firdes_kaiser_window interface.

#include <liquid/liquid.h>

int main() {
    // options
    float fc=0.15f;         // filter cutoff frequency
    float ft=0.05f;         // filter transition
    float As=60.0f;         // stop-band attenuation [dB]
    float mu=0.0f;          // fractional timing offset

    // estimate required filter length and generate filter
    unsigned int h_len = estimate_req_filter_len(ft,As);
    float h[h_len];
    liquid_firdes_kaiser(h_len,fc,As,mu,h);
}

An example of a low-pass filter design using the Kaiser window can be found in [ref:fig-filter-firdes_kaiser] .

Nyquist filter design∞

Nyquist's criteria for designing a band-limited filter without inter-symbol interference is for the spectral response of a linear phase filter to be symmetric about its symbol rate. liquid provides several Nyquist filters as design prototypes and are listed in [ref:tab-filter-nyquist] .

Table [tab-filter-nyquist]. Nyquist filter prototypes available in liquid

Scheme	Description
LIQUID_NYQUIST_KAISER	Kaiser filter
LIQUID_NYQUIST_PM	Parks-McClellan algorithm
LIQUID_NYQUIST_RCOS	Raised-cosine filter
LIQUID_NYQUIST_FEXP	flipped exponential {cite:Beaulieu:2001}
LIQUID_NYQUIST_FSECH	flipped hyperbolic secant {cite:Assalini:2004}
LIQUID_NYQUIST_FARCSECH	flipped hyperbolic arc-secant {cite:Assalini:2004}

The interface for designing Nyquist filters is simply

liquid_firdes_nyquist(_ftype, _k, _m, _beta, _dt, *h);

where _ftype is one of the filter types in[ref:tab-filter-nyquist] ,\(k\) is the number of samples per symbol,\(m\) is the filter delay in symbols,\(\beta\) is the excess bandwidth (rolloff) factor,\(\Delta t\) is the fractional sample delay (usually set to zero for typical filter designs and is ignored in the LIQUID_NYQUIST_PM design), and \(h\) is the output coefficients array of length \(2km+1\) .

Square-Root Nyquist filter design∞

Square-root Nyquist filters are commonly used in digital communications systems with linear modulation as a pulse shape for matched filtering. Applying a pulse shape to the transmitted symbol sequence limits its occupied spectral bandwidth by smoothing the transitions between symbols. If an identical filter is applied at the receiver then the system is matched resulting in the maximum signal-to-noise ratio and (theoretically) zero inter-symbol interference. While the design of Nyquist filters is trivial and can be accomplished by applying any desired window to a \(\sinc\) function, designing a square-root Nyquist filter is not as straightforward. liquid conveniently provides several square-root Nyquist filter prototypes listed in [ref:tab-filter-rnyquist] .

Table [tab-filter-rnyquist]. Square-root Nyquist filter prototypes available in liquid

Scheme	Description
LIQUID_RNYQUIST_ARKAISER	Approximate square-root Kaiser filter
LIQUID_RNYQUIST_RKAISER	Square-root Kaiser
LIQUID_RNYQUIST_RRCOS	Square-root raised-cosine filter
LIQUID_RNYQUIST_hM3	harris-Moerder type 3 {cite:harris-Moerder:2005}
LIQUID_RNYQUIST_GMSTX	GMSK transmit filter {cite:Proakis:2001}
LIQUID_RNYQUIST_GMSRX	GMSK receive filter
LIQUID_RNYQUIST_FEXP	flipped exponential {cite:Beaulieu:2001}
LIQUID_RNYQUIST_FSECH	flipped hyperbolic secant {cite:Assalini:2004}
LIQUID_RNYQUIST_FARCSECH	flipped hyperbolic arc-secant {cite:Assalini:2004}

The interface for designing square-root Nyquist filters is simply

liquid_firdes_rnyquist(_ftype, _k, _m, _beta, _dt, *h);

where _ftype is one of the filter types in [ref:tab-filter-rnyquist] ,\(k\) is the number of samples per symbol,\(m\) is the filter delay in symbols,\(\beta\) is the excess bandwidth (rolloff) factor,\(\Delta t\) is the fractional sample delay (usually set to zero for typical filter designs), and \(h\) is the output coefficients array of length \(2km+1\) . All square-root Nyquist filters in liquid have these four basic properties (\(k\) , \(m\) , \(\beta\) , \(\Delta t\) ) and produce a filter with\(N=2km+1\) coefficients.

The most common square-root Nyquist filter design in digital communications is the square-root raised-cosine (RRC) filter, likely due to the fact that an expression for its time series can be expressed in closed form. The filter coefficients themselves are derived from the following equation:

$$ h\left[z\right] = 4\beta \frac{ \cos\left[(1+\beta)\pi z\right] + \sin\left[(1-\beta)\pi z\right] / (4\beta z) } { \pi \sqrt{T}\left[ 1-16\beta^2z^2\right] } $$

where \(z=n/k-m\) , and \(T=1\) for most cases. liquid compensates for the two cases where \(h[n]\) might be undefined in the above equation, i.e.

$$ \mathop {\lim }\limits_{z \to 0 } h(z) = 1 - \beta + 4\beta/\pi $$

and

$$ \mathop {\lim }\limits_{z \to \pm \frac{1}{4\beta} } h(z) = \frac{\beta}{\sqrt{2}} \left[ \left(1 + \frac{2}{\pi}\right)\sin\left(\frac{\pi}{4\beta}\right) + \left(1 - \frac{2}{\pi}\right)\cos\left(\frac{\pi}{4\beta}\right) \right] $$

The \(r\) -Kaiser and harris-Moerder-3 (hM3) filters cannot be expressed in closed form but rely on iterations over traditional filter design techniques to search for the design parameters which minimize the resulting filter's inter-symbol interference (ISI). Similarly the approximate \(r\) -Kaiser filter uses an approximation for the design parameters to eliminate the need for running the search; this comes at the expense of a slight performance degradation.

[ref:fig-filter-rnyquist] contrasts the different square-root Nyquist filters available in liquid. The square-root raised-cosine filter is inferior to the (approximate) \(r\) -Kaiser and harris-Moerder-3 filters in both transition bandwidth as well as side-lobe suppression. In the figure the responses of the \(r\) -Kaiser and approximate \(r\) -Kaiser filters are indistinguishable.

GMSK Filter Design∞

The transmit filter for a GMSK modem with a bandwidth-time product \(BT\) (equivalent to the excess bandwidth factor \(\beta\) ) is defined as

$$ h_t(t) = Q\left( \frac{2 \pi BT}{\sqrt{\ln(2)}} \left(t-\frac{1}{2}\right) \right) - Q\left( \frac{2 \pi BT}{\sqrt{\ln(2)}} \left(t+\frac{1}{2}\right) \right) $$

where \(Q(z)=\frac{1}{2}\left(1 - erf(z/\sqrt{2})\right)\) (see [ref:section-math-transcendentals-Q] ). The transmit filter imparts inter-symbol interference, leaving the receiver to compensate. liquid implements a GMSK receive filter by minimizing the inter-symbol interference of the composite, and as such there is no closed-form solution for the GMSK receive filter.

[ref:fig-filter-firdes_gmskrx] depicts the transmit, receive, and composite filters in both the time and frequency domains. Notice that the frequency response of the receive filter has a gain in the pass-band (around \(f=0.13\) ) to compensate for the ISI imparted by the transmit filter. Consequently the composite filter has nearly zero ISI, as can be seen by its flat pass-band response and transition through \(20\log_{10}\left(\frac{1}{2}\right)\) .

Parks-McClellan Algorithm∞

FIR filter design using the Parks-McClellan algorithm is implemented in liquid with the firdespm interface. The Parks-McClellan algorithm uses the Remez exchange algorithm to solve the minimax problem (minimize the maximum error) for filter design. The interface accepts a description of \(N_b\) disjoint and non-overlapping frequency bands with a desired response and relative error weighting for each, and computes the resulting filter coefficients.

firdespm_run(_h_len,        // filter length
             _bands*,       // array of frequency bands
             _des*,         // desired response in each band
             _weights*,     // relative weighting for each band
             _num_bands,    // number of bands
             _btype,        // filter type
             _wtype*,       // weighting function for each band
             _h*)

• _bands is a \([N_b \times 2]\) matrix of the band edge descriptions. Each row corresponds to an upper and lower band edge for each region of interest. These regions cannot be overlapping.
• _des is an array of size \(N_b\) with the desired response (linear) for each band.
• _weights is an array of size \(N_b\) with the relative error weighting for each band.
• _num_bands represents \(N_b\) , the number of bands in the design.
• _btype gives the filter type for the design. This is typically LIQUID_FIRDESPM_BANDPASS for the majority of filters.
• _wtype is an array of length \(N_b\) which specifies the weighting function for each band (flat, exponential, or linear).

Listed below is an example of the firdespm interface.

#include <liquid/liquid.h>

int main() {
    // define filter length, type, number of bands
    unsigned int n=55;
    liquid_firdespm_btype btype = LIQUID_FIRDESPM_BANDPASS;
    unsigned int num_bands = 4;

    // band edge description [size: num_bands x 2]
    float bands[8]   = {0.0f,   0.1f,   // 1    first pass-band
                        0.15f,  0.3f,   // 0    first stop-band
                        0.33f,  0.4f,   // 0.1  second pass-band
                        0.42f,  0.5f};  // 0    second stop-band

    // desired response [size: num_bands x 1]
    float des[4]     = {1.0f, 0.0f, 0.1f, 0.0f};

    // relative weights [size: num_bands x 1]
    float weights[4] = {1.0f, 1.0f, 1.0f, 1.0f};

    // in-band weighting functions [size: num_bands x 1]
    liquid_firdespm_wtype wtype[4] = {LIQUID_FIRDESPM_FLATWEIGHT,
                                      LIQUID_FIRDESPM_EXPWEIGHT,
                                      LIQUID_FIRDESPM_EXPWEIGHT,
                                      LIQUID_FIRDESPM_EXPWEIGHT};

    // allocate memory for array and design filter
    float h[n];
    firdespm_run(n,num_bands,bands,des,weights,wtype,btype,h);
}

Miscellaneous functions∞

Here are several miscellaneous functions used in liquid's filter module, useful to filtering and filter design.

• estimate_req_filter_len(df,As) returns an estimate of the required filter length, given a transition bandwidth \(\Delta f\) and stop-band attenuation \(A_s\) . The estimate uses Kaiser's formula {cite:Vaidyanathan:1993}

$$ N \approx \frac{ A_s - 7.95 }{ 14.26 \Delta f } $$

• estimate_req_filter_As(df,N) returns an estimate of the filter's stop-band attenuation \(A_s\) given the filter's length \(N\) and transition bandwidth \(\Delta f\) . The estimate uses an iterative binary search to find \(A_s\) from estimate_req_filter_As() .

• estimate_req_filter_df(As,N) returns an estimate of the filter's transition bandwidth \(\Delta f\) given the filter's length \(N\) and stop-band attenuation \(A_s\) . The estimate uses an iterative binary search to find \(\Delta f\) from estimate_req_filter_As() .

• kaiser_beta_As(As) returns an estimate of the Kaiser \(\beta\) factor for a particular stop-band attenuation \(A_s\) . The estimate uses Kaiser's original formula {cite:Vaidyanathan:1993}, viz

$$ \beta = \begin{cases} 0.1102 (A_s - 8.7) & A_s \gt 50 \\ 0.5842 (A_s - 21)^{0.4} & 21 \lt A_s \le 50 \\ 0 & \text{else} \end{cases} $$

• fir_group_delay(*h,n,f) computes the group delay for a finite impulse-response filter with\(n\) coefficients \(\vec{h}\) at a frequency \(f\) . The group delay \(\tau_g\) at frequency \(f\) for a finite impulse response filter of length \(N\) is computed as

$$ \tau_g = \Re\Biggl\{ \,\, \frac{ \sum_{k=0}^{N-1}{h(k)e^{j 2 \pi f k} \cdot k} } { \sum_{k=0}^{N-1}{h(k)e^{j 2 \pi f k}} }\,\, \Biggr\} $$

• iir_group_delay(*b,nb,*a,na,f) computes the group delay for an infinite impulse-response filter with \(n_a\) feed-back coefficients \(\vec{a}\) , and \(n_b\) feed-forward coefficients \(\vec{b}\) at a frequency \(f\) . To compute the group delay, define the flipped convolution of \(\vec{a}\) and \(\vec{b}\) as\(c(n) = \sum_{m=0}^{N-1}{a^*(m) b(m-n)}\) for \(n \in \{0,1,\ldots,2(N+1)\}\) . The group delay \(\tau_g\) at frequency \(f\) for an infinite impulse response filter of order \(N\) can therefore be computed as

$$ \tau_g = \Re\Biggl\{ \,\, \frac{ \sum_{k=0}^{2(N+1)}{c(k)e^{j 2 \pi f k} \cdot k} } { \sum_{k=0}^{2(N+1)}{c(k)e^{j 2 \pi f k}} }\,\, \Biggr\} - N $$

• iirdes_isstable(*b,*a,n) checks the stability of an infinite impulse-response filter with \(n\) feed-back and feed-forward coefficients \(\vec{a}\) and \(\vec{b}\) respectively. Stability is tested by computing the roots of the denominator (poles) and ensuring that they lie within the unit circle. Notice that the poles in[ref:fig-filter-butter] -[ref:fig-filter-bessel] all have their poles within the unit circle and are therefore stable (as expected).

• liquid_filter_autocorr(*h,N,n) computes the auto-correlation of a filter with an array of coefficients \(\vec{h}\) of length \(N\) at a specific lag \(n\) as

$$ r_{hh}(n) = \sum_{k=n}^{N-1} {h(k)h^*(k-n)} $$

• liquid_filter_isi(*h,k,m,*rms,*max) computes the inter-symbol interference (both mean-squared error and maximum error) for a filter \(\vec{h}\) with \(k\) samples per symbol and delay of \(m\) samples. The filter has \(2km+1\) coefficients and the resulting RMS and maximum ISI are stored in rms and max , respectively. This is useful in comparing the performance of root-Nyquist matched filter designs (e.g. root raised-cosine).

• liquid_filter_energy(*h,N,fc,nfft) computes the relative out-of-band energy \(E_0\) at a cutoff frequency\(f_c\) for a finite impulse response filter \(\vec{h}\) with \(N\) coefficients. The parameter nfft specifies the precision of the computation. The relative out-of-band energy is computed as

$$ E_0 = \frac{ \int_{2 \pi f_c}^{\infty}{|H(\omega)|^2d\omega} } { \int_{0}^{\infty}{|H(\omega)|^2d\omega} } $$