StackOverflow2013

Note that there are some explanatory texts on larger screens.

plurals

PO
text
Body
copied!<p>Tom's answer is the following recursive filter:</p> <pre><code>y[n] = (1 - a)*y[n-1] + a*x[n] H(z) = Y(z)/X(z) = a / (1 - (1 - a)*1/z) </code></pre> <p>I'll plot this in Python/pylab for a=0.25, a=0.50, and a=0.75:</p> <pre><code>from pylab import * def H(a, z): return a / (1 - (1 - a) / z) w = r_[0:1000]*pi/1000 z = exp(1j*w) H1 = H(0.25, z) H2 = H(0.50, z) H3 = H(0.75, z) plot(w, abs(H1), 'r') # red plot(w, abs(H2), 'g') # green plot(w, abs(H3), 'b') # blue </code></pre> <p><img src="https://i.stack.imgur.com/3pRvS.png" alt="alt text"></p> <p>Pi radians/sample is the Nyquist frequency, which is half the sampling frequency. </p> <p>If this simple filter is inadequate, try a 2nd order Butterworth filter:</p> <pre><code># 2nd order filter: # y[n] = -a[1]*y[n-1] - a[2]*y[n-2] + b[0]*x[n] + b[1]*x[n-1] + b[2]*x[n-2] import scipy.signal as signal # 2nd order Butterworth filter coefficients b,a # 3dB cutoff = 2000 Hz fc = 2000.0/44100 b, a = signal.butter(2, 2*fc) # b = [ 0.01681915, 0.0336383 , 0.01681915] # a = [ 1. , -1.60109239, 0.66836899] # approximately: # y[n] = 1.60109*y[n-1] - 0.66837*y[n-2] + # 0.01682*x[n] + 0.03364*x[n-1] + 0.01682*x[n-2] # transfer function def H(b,a,z): num = b[0] + b[1]/z + b[2]/(z**2) den = a[0] + a[1]/z + a[2]/(z**2) return num/den H4 = H(b, a, z) plot(w, abs(H4)) # show the corner frequency plot(2*pi*fc, sqrt(2)/2, 'ro') xlabel('radians') </code></pre> <p><img src="https://i.stack.imgur.com/IdWVn.png" alt="alt text"></p> <p>Evaluate a test signal at the 3dB cutoff frequency <code>fc=2000</code>:</p> <pre><code>fc = 2000.0/44100 b, a = signal.butter(2, 2*fc) # test signal at corner frequency (signed 16-bit) N = int(5/fc) # sample for 5 cycles x = int16(32767 * cos(2*pi*fc*r_[0:N])) # signed 16-bit output yout = zeros(size(x), dtype=int16) # temp floats y = 0.0 y1 = 0.0 y2 = 0.0 # filter the input for n in r_[0:N]: y = (-a[1] * y1 + -a[2] * y2 + b[0] * x[n] + b[1] * x[n-1] + b[2] * x[n-2]) # convert to int16 and saturate if y > 32767.0: yout[n] = 32767 elif y < -32768.0: yout[n] = -32768 else: yout[n] = int16(y) # shift the variables y2 = y1 y1 = y # plots plot(x,'r') # input in red plot(yout,'g') # output in green # show that this is the 3dB point plot(sqrt(2)/2 * 32768 * ones(N),'b-') xlabel('samples') </code></pre> <p><img src="https://i.stack.imgur.com/7O1qX.png" alt="alt text"></p>

Querying!

Guidance

An individual column

Larger individual text columns get their own page to allow for proper reading.

SQuiL has stopped working due to an internal error.

If you are curious you may find further information in the browser console, which is accessible through the devtools (F12).

Reload