3

I have created an IIR design algorithm that generates complex coefficients (there is no symmetry in the poles and zeros). However, the IIRs will be used to filter a real signal. Is there a closed form (analytical) way to convert the IIR to real coefficients. I did think that using the output of the difference equation to solve a system of linear equations might be a known solution?

Another way to put the question: when both complex and real IIRs are acting as a digital filter processing discrete samples, for a given real input they produce the same real output (discarding the imaginary component).

Another way to put the question: Both the complex and real IIRs have the same magnitude response.

keith
  • 906
  • 1
  • 7
  • 20
  • converting a real signal to complex would seem more straightforward –  Sep 20 '17 at 13:03
  • @Stanley Pawlukiewicz, agreed. I'm looking for insight that would reduce the number of operations per sample to process the filter i.e. if an 8th order complex IIR can easily be transformed to an 8th order real IIR, the operation count drops (and the state storage halves). – keith Sep 20 '17 at 13:26
  • I don't see how you can get a real polynomial without conjugate symmetric roots –  Sep 20 '17 at 13:33
  • @Stanley Pawlukiewicz, I suspect you can't help me. The question is an analytical way to convert a complex IIR to a real IIR. An obvious a side effect of such a conversion would be a rational polynomial with a mixture of real and conjugate pairs of poles and zeros. – keith Sep 20 '17 at 13:38
  • @Olli Niemitalo, for a given real input it produces the same real output when processing discrete samples when acting as a digital filter. – keith Sep 20 '17 at 13:39
  • 4
    The real output is half the sum of the complex output and its complex conjugate. The complex conjugate output is produced by the filter with complex conjugate coefficients. Using linearity, the filter you want is the parallel concatenation of the original filter and the complex conjugate filter, scaled by 1/2. – Jazzmaniac Sep 20 '17 at 13:44
  • @Jazzmaniac, that makes sense, so simple, thanks :-) Shame the order doubles - I guess there's no free lunch and the operation count stays the same. – keith Sep 20 '17 at 13:49
  • 1
    @Jazzmaniac what is parallel concatenation? – Olli Niemitalo Sep 20 '17 at 13:59
  • @Olli Niemitalo, he means sum the finite difference output and divide by two. That's the easiest way to look at it anyway. Maybe he'll elaborate in an real answer and I can mark it correct :-) – keith Sep 20 '17 at 14:01
  • 1
    Your new edits are problematic. It is important to mention that you discard the imaginary part of the filter output. – Olli Niemitalo Sep 22 '17 at 11:20
  • 1
    @Olli Niemitalo, yes I feel embarrassed that I haven't defined the problem particularly well. I am not sure if I need to defined the output of the complex filter. I'm hoping the impulse response of the real filter is well defined by the change in constraint? If it doesn't make sense, I probably need some help understanding why not. – keith Sep 22 '17 at 11:53
  • Non-zero imaginary parts in any of the coefficients of the complex IIR filter will break the Hermitian symmetry of the frequency response (See Fig. 1 in my answer). Frequency responses $F(\omega)$ of all IIR filters with real coefficients have Hermitian symmetry $F(\omega) = F^*(-\omega).$ Hermitian symmetry is manifested in the magnitude frequency response of real-coefficient filters as $|F(\omega)| = |F(-\omega)|.$This makes the last sentence of the question a contradiction, in the general sense. – Olli Niemitalo Sep 22 '17 at 17:50

1 Answers1

5

This answer shows how to create a real-coefficient infinite-impulse-response (IIR) filter, the output of which equals the real part of the output of a given complex-coefficient IIR filter. Also an example is given that shows that the approach can be useful.

Z-transform $\mathcal{Z}\{h[n]\}$ of a time-domain sequence $h[n]$ has the property that:

$$\mathcal{Z}\{h^*[n]\} = H^*(z^*),\quad\text{where }H(z) = \mathcal{Z}\{h[n]\}.\tag{1}$$

Given an IIR filter transfer function:

$$H(z) = H_0\frac{\prod^M_{k=1}\left(1-q_kz^{-1}\right)}{\prod^N_{k=1}\left(1-p_kz^{-1}\right)},\tag{2}$$

with real coefficient $H_0$ and $M$ complex zeros $q_k$ and $N$ complex poles $p_k,$ discarding of the imaginary part of its output $h[n]$ to form real output $g[n]$:

$$g[n] = \operatorname{Re}(h[n]) = \frac{1}{2}\left(h[n]+h^*[n]\right),\tag{3}$$

is equivalent to having a filter with transfer function:

$$G(z) = \frac{1}{2}\big(H(z) + H^*(z^*)\big).\tag{4}$$

Equivalently, with $z = e^{i\omega},$ this new filter has a frequency response $\frac{1}{2}\big(F(\omega) + F^*(-\omega)\big),$ where $F(\omega)$ is the frequency response of the full complex filter. Using further properties of complex conjugation:

$$G(z) = \frac{1}{2}H_0\left(\frac{\prod^M_{k=1}\left(1-q_kz^{-1}\right)}{\prod^N_{k=1}\left(1-p_kz^{-1}\right)} + \frac{\prod^M_{k=1}\left(1-q^*_kz^{-1}\right)}{\prod^N_{k=1}\left(1-p^*_kz^{-1}\right)}\right)\\ = \frac{\frac{1}{2}H_0\left(\prod^M_{k=1}\left(1-q_kz^{-1}\right)\prod^N_{k=1}\left(1-p^*_kz^{-1}\right) + \prod^M_{k=1}\left(1-q^*_kz^{-1}\right)\prod^N_{k=1}\left(1-p_kz^{-1}\right)\right)}{\prod^N_{k=1}\left(1-p_kz^{-1}\right)\prod^N_{k=1}\left(1-p^*_kz^{-1}\right)}\\ = \frac{\mathcal{Z}\bigg\{\operatorname{Re}\Big(\mathcal{Z}^{-1}\left\{H_0\prod^M_{k=1}\left(1-q_kz^{-1}\right)\prod^N_{k=1}\left(1-p^*_kz^{-1}\right)\right\}\Big)\bigg\}}{\prod^N_{k=1}\Big(\left(1-p_kz^{-1}\right)\left(1-p^*_kz^{-1}\right)\Big)}.$$$\tag{5}$

By expanding the products to polynomials of $z^{-1}$, the above gives the procedure to calculate the coefficients of the real filter given the coefficients of the original complex filter, by convolutions of sequences of polynomial coefficients. If the original complex filter with input $x[n]$ and output $y[n]$ is in form:

$$\begin{align}y[n] &= b_0\,x[n] + b_1\,x[n-1] + b_2\,x[n-2] + \ldots\\ &-a_1\,y[n-1] - a_2\,y[n-1] - \ldots,\end{align}\tag{6}$$

that readily accepts the coefficients $B = [b_0,\,b_1,\,\ldots]$ and $A = [a_0,\,a_1,\,\ldots]$ of the numerator and denominator polynomials $b_0+b_1\,z^{-1}+b_2\,z^{-2}+\ldots$ and $a_0+a_1\,z^{-1}+a_2\,z^{-2}+\ldots,$ respectively, with a hidden coefficient $a_0 = 1$, then the coefficients $B_\text{real}$ and $A_\text{real}$ for the new filter will be given by:

$$\begin{align}B_\text{real} &= \frac{1}{2}\left(B*A^* + B^**A\right) = \operatorname{Re}\left( B*A^*\right)\\ A_\text{real} &= A*A^* = \operatorname{Re}\left(A*A^*\right),\end{align}\tag{7}$$

where $*$ represents convolution (or complex conjugate if superscript). The last $\operatorname{Re}$ is unnecessary, but may be useful in numerical convolution.

Example 1: Minimal example

A minimal example in Octave starts with a complex filter and generates a real filter that matches the real output of the complex filter:

b = [0.5, 0.35 - 0.35*j]; #Coefficients of numerator polynomial
a = [1, - 0.6 - 0.6*j]; #Coefficients of denominator polynomial
length = 1024;
impulse = eye(1, length);
[h_complex, w_complex] = freqz(b, a, length, "whole");
[h_discard_imag, w_discard_imag] = freqz(real(filter(b, a, impulse)), [1], length, "whole");
b_real = real(conv(b, conj(a)))
a_real = real(conv(a, conj(a)))
[h_real, w_real] = freqz(filter(b_real, a_real, impulse), [1], length, "whole");
freqz_plot([w_complex, w_discard_imag, w_real], [h_complex, h_discard_imag, h_real]);

The generated real coefficients of the numerator and denominator polynomials of the new real filter are:

b_real =
   0.50000   0.05000   0.00000

a_real =
   1.00000  -1.20000   0.72000

The frequency response of the new real filter equals that calculated from the real part of the impulse response of the complex filter:

Frequency responses
Figure 1. Blue: Frequency response of the complex filter with the imaginary part of its output still present. The frequency response lacks Hermitian symmetry, so no real-coefficient filter can reproduce it. Red: Shared frequency response of the new real-coefficient filter and the complex-coefficient filter with the imaginary part of its output discarded. The shared frequency response has Hermitian symmetry.

Example 2: Sliding discrete Fourier transform

Constructing a real filter from a complex filter can sometimes be useful. One such case is calculation of the real part of a single term of a sliding discrete Fourier transform (sliding DFT) of length $r$ at frequency $\omega = 2\pi m/r,\, m\in N.$ The complex IIR filter is:

$$\begin{align}y[n] &= x[n] - x[n-r]\\ &+e^{-i\omega}y[n-1]\end{align}\tag{8}$$

We get using Eqs. 6 and 7 the coefficients of the numerator and denominator polynomials of the transfer function of the real filter $(A_\text{real},\,B_\text{real})$ from those of the complex filter $(A,\,B).$ Assuming that $r > 1$ (coefficients not given are zero):

$$\begin{array}{rcllll} k&=&0&1&r&r+1\\ \hline B[k]&=&1&&-1\\ A[k]&=&1&-e^{-i\omega}\\ B_\text{real}[k]&=&\operatorname{Re}(1)&\operatorname{Re}(-e^{i\omega})&\operatorname{Re}(-1)&\operatorname{Re}(e^{i\omega})\\ &=&1&-\cos(\omega)&-1&\cos(\omega)\\ A_\text{real}[k]&=&1&-e^{i\omega}-e^{-i\omega}&-e^{-i\omega}(-e^{i\omega})\\ &=&1&-2\cos(\omega)&1\end{array}$$

The real filter is thus:

$$\begin{align}y[n] &= x[n] -\cos(\omega)x[n-1] -x[n-r] -\cos(\omega)x[n-r-1]\\ &+ 2\cos(\omega)y[n-1] - y[n-2]\end{align}$$

The recursive part can be recognized as the Goertzel filter.

In Octave:

m = 2;
r = 20;
w = 2*pi*m/r;
b = zeros(1, r+2);
b(1) = 1; b(2) = -cos(w); b(r+1) = -1; b(r+2) = cos(w);
a = [1, -2*cos(w), 1];
length = r+10;
impulse = eye(1, length);
plot([0:length-1], filter(b, a, impulse), "o");

Real part of impulse response of the real filter
Figure 2. Impulse response of a real IIR filter that calculates the real part of a single frequency bin of sliding DFT.

The impulse response can be zero phase aligned to the opposite end by multiplying the complex transfer function by $e^{-i\omega}:$

m = 2;
r = 20;
w = 2*pi*m/r;
b = zeros(1, r+2);
b(1) = cos(w); b(2) = -1; b(r+1) = -cos(w); b(r+2) = 1;
a = [1, -2*cos(w), 1];
length = r+10;
impulse = eye(1, length);
plot([0:length-1], filter(b, a, impulse), "o");

Impulse response with opposite zero phase alignment
Figure 3. Impulse response of a real IIR filter that calculates the real part of a single frequency bin of sliding DFT, with zero phase aligned to the opposite end than in Fig. 2. This is a circular shift of the impulse response one sampling period to the left, and can sometimes be used to remove one sample of system latency.

In case of $m=0,$ the numerator and denominator polynomials will have an unnecessary common factor, but for other values of $m$ the approach works fine.

Olli Niemitalo
  • 13,491
  • 1
  • 33
  • 61
  • I need to double check this before I mark it correct, but much appreciate the time you've taken to explain it so well. – keith Sep 22 '17 at 10:16
  • I'm going to check with some designs that contain single complex zeros on the unit disk (such that the angular frequency of the zero is negative) - empirically I've found this to be troublesome. – keith Sep 22 '17 at 10:22
  • The magnitude response of $G(z)$ is not equal to the magnitude response of $H(z)$, which means the approach produces real output, but not of the same magnitude response as the original filter. I think for the approach to work the condition is actually $H^(z) = H(z^)$ which implies conjugate symmetry of the response and therefore of the complex poles/zeros and therefore produces a real coefficient filter anyway. Then a $G(z)$ can be constructed with the same magnitude response as $H(z)$. I upvoted you anyways :-) – keith Sep 22 '17 at 11:02
  • Try these coefficients: a = [1, -0.2266324680203 + 0.973980351156536*j]; and b = [1]; – keith Sep 22 '17 at 11:58
  • @keith Those a and b also give equal plots for _discard_imag and _real. The notch in the magnitude frequency response of the _complex filter (also present in my example) disappears in _discard_imag because discarding the imaginary part gives a new frequency response $\frac{1}{2}(F(\omega) + F^*(-\omega))$ from the original frequency response $F(\omega)$ and the notch does not survive addition of large-magnitude stuff to it. – Olli Niemitalo Sep 22 '17 at 12:16
  • Olli the answer nicely illustrates the symmetry properties of LTI systems and it's correct however the last line would be $B_{real} = 0.5 ( b \star a^{} + a \star b^{} ) $ and $A_{real} = a \star a^{*}$, @keith , give it a try. But of course it does not reduce the number of arithmetics as the order is doubled. However I still doubt if there could be another method to really reduce the count by two, as in the first place half of the computations is thrown away... – Fat32 Sep 22 '17 at 19:53
  • My $B$ and $A$ hold place for the general $a$ (denominator) and $b$ (numerator) coefficients, I thought yours was also so. In addition, I assumed the final filter to be $G(z)$ which produces the real part of the complex output that's produced by the original complex filter $a$ and $b$ (which is produced by the real part of the impulse response $h[n]$) – Fat32 Sep 22 '17 at 20:10
  • Yes the edit now clarifies they are the same. I wonder why you preferred the reverse... It could be better if you have written explicitly as numerator and denominator coefficients next to $A$ and $B$ respectively for those who read the site while multitasking with other things :-) – Fat32 Sep 22 '17 at 20:19
  • for what (reasonable) reason ? – Fat32 Sep 22 '17 at 20:30
  • hmm plausible. The left side of the LCCDE is also read first so I use $a$ there :-). – Fat32 Sep 22 '17 at 20:32
  • Octave, Matlab and (many) DTSP books use the opposite (to yours) notation. It's just a notation of course. I admit while reading your answer I skipped the difference equation part which indicated your choice explicitly. Anyway, it's quite interesting that you have learned Octave/Matlab notation today. So you don't use Matlab/Octave since 98? – Fat32 Sep 22 '17 at 20:48
  • Yes they refer to the same thing here but I think divident and divisor refer to the actual operation of division (which also has the quotient and remainder outputs) while the numerator, denominator refers to the fraction symbolism. – Fat32 Sep 22 '17 at 21:07
  • I have now swapped a and b in the notation to reflect the more common use. – Olli Niemitalo Sep 23 '17 at 06:13