0
$\begingroup$

I've been trying to design a bandpass filter using scipy but I keep getting a LinAlg Singular Matrix error. I read that a singular matrix is one that is not invertable, but I'm not sure how that error is coming up and what I can do to fix it

The code takes in an EEG signal (which, in the code below, I have just replaced with an int array for testing) and filters out frequencies < 8Hz and > 12Hz (alpha band)

Can anyone shed some light on where the singular matrix error is coming from? Or alternatively, if you know of a better way to filter a signal like this I'd love to test out other options too

from scipy import signal
from scipy.signal import filter_design as fd
import matplotlib.pylab as plt

#bandpass
Wp = [8, 12]   # Cutoff frequency
Ws = [7.5, 12.5]   # Stop frequency
Rp = 1             # passband maximum loss (gpass)
As = 100              # stoppand min attenuation (gstop)

b,a = fd.iirdesign(Wp,Ws,Rp,As,ftype='butter')
w,H = signal.freqz(b,a)  # filter response
plt.plot(w,H)

t = np.linspace(1,256,256)
x = np.arange(256)
plt.plot(t,x)

y = signal.filtfilt(b,a,x)
plt.plot(t,y)
Improve this question
$\endgroup$
2

2 Answers 2

Reset to default
1
$\begingroup$

Sometimes Stackexchange coughs up these ancient unanswered questions. Since I've started, I'll answer in spite of the age:

iirdesign expects you to give it frequencies in the range 0..1, normalized to the Nyquist frequency (the latest version looks like it'll let you tell it the sampling rate for your own convenience). You can expect that any filter design software that's designing a sampled-time frequency and isn't given the sampling rate is going to expect the frequency parameters to be normalized somehow -- as iirdesign does, or normalized to the sampling rate, or as radians/sample, etc.

If you're using unfamiliar design software, you should (A) be aware of this and read the documentation to see how you're supposed to be scaling things, and (B) read the documentation anyway, because there's always unexpected details (for instance, I would have expected frequencies to be in radians/sample, not proportions of Nyquist, so I would have given it the wrong numbers).

Probably what happened here is that iirdesign got nonsensical frequencies and either coughed up the singular matrix error, or it did its best, came up with a filter, and filtfilt coughed up the error.

Improve this answer
$\endgroup$
0
$\begingroup$

I'd suggest checking again your parameter choices for

b,a = fd.iirdesign(Wp,Ws,Rp,As,ftype='butter')

Looking at the resulting filter coefficients, you have pole/zero cancellation at DC, i.e. at $z=1$.

Improve this answer
$\endgroup$
4
  • $\begingroup$ it turns out that iirdesign requires both Wp and Ws to be normalised. That ended up fixing my problem. Having said that though, I'm very new to signal processing - could you explain what you mean by pole cancellation at DC? $\endgroup$
    – Simon
    Commented Jul 31, 2015 at 20:02
  • $\begingroup$ The $z$-transform of the impulse response produced by your code has the form $$H(z) = K \frac{(1-az^{-1})(1-bz^{-1})}{(1-a^{-1})(1-c^{-1})}$$ where $K$ is the system gain. The parametersin the numerator are referred to as zeros, e.g. $H(z) = 0$ for $z=b$. Likewise, the parameters in the denominator are referred to as poles since they make $H(z)$ blow up, i.e.$H(c) = \infty$. So, with $a$ being both a pole and zero we have that $H(a) = \frac{0}{0}$ which is only an issue numerically since an equivalent $H(z)$ is $$H(z) = K\frac{(1-bz^{-1})}{(1-cz^{-1})}$$ which does not have this problem. $\endgroup$
    – curiousStudent
    Commented Jul 31, 2015 at 22:21
  • $\begingroup$ hmm so am I right in thinking that cancellation is a bad thing? I've tried to read up on this but it is far beyond my level of understanding. Are we trying to prevent a case where both poles and zeros = 1? $\endgroup$
    – Simon
    Commented Aug 2, 2015 at 0:28
  • 1
    $\begingroup$ It can be, it all depends on the implementation details. Essentially it comes down to the difference between IO (input output) stability and internal stability. Thnk about computing $\frac{x^2}{x}$ for $x=0$. One the one hand, we have $\frac{0}{0}$ from which we can't conclude anything. On the other hand, we know the answer is $0$ since $\frac{x^2}{x}=x$. So, a naive implementation would be to compute the numerator, compute the denominator, and then divide and would result in an error. Checking for cancellation first eliminates this problem and would give the correct answer. $\endgroup$
    – curiousStudent
    Commented Aug 2, 2015 at 17:24

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.