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)

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$.

  • 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? – Simon Jul 31 '15 at 20:02
  • 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. – curiousStudent Jul 31 '15 at 22:21
  • 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? – Simon Aug 2 '15 at 0:28
  • 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. – curiousStudent Aug 2 '15 at 17:24

Your Answer

 
discard

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.