In [2]:
import scipy.signal

In [8]:
def plotspec(x, Ts):
fig = figure()
ax1.plot(x)

q = fft.fft(x)
ax2.plot(fft.fftfreq(len(x), Ts), abs(q))

time = 3.0
Ts=1.0/10000.0
x = random.uniform(-1.0,1.0, time/Ts)
plotspec(x, Ts)

figure()
b = scipy.signal.remez(100, [0,0.1,0.105,0.5],[1,0])
y = scipy.signal.lfilter(b,1,x)
plotspec(y, Ts)

figure()
b = scipy.signal.remez(100, [0,0.12,0.13,0.25,0.255, 0.5],[0,1,0])
y = scipy.signal.lfilter(b,1,x)
plotspec(y, Ts)

figure()
b = scipy.signal.remez(100, [0,0.37,0.38,0.5],[0,1])
y = scipy.signal.lfilter(b,1,x)
plotspec(y, Ts)


Notes that remez() takes slightly different arguments in scipy.signal than it does in MATLAB:

• Top end of frequency specification is 0.5 instead of 1.0 (ie: Nyquist frequency), so all those terms are halved
• The third argument takes half as many inputs as the previous one: it expects that we're defining levels between the bands, or something

3.8. Mimic the code in filternoise to create a filter that

1. passes all frequencies above 500Hz
2. passes all frequencies below 3000Hz
3. rejects all frequences between 1500 and 2500 Hz
In [20]:
def plotfilter(b, Ts):
figure()
time = 3.0
x = random.uniform(-1.0,1.0, time/Ts)
y = scipy.signal.lfilter(b,1,x)
plotspec(y, Ts)

def genlpf(ntaps, topf, Ts):
q = topf/(1/Ts)
b = scipy.signal.remez(ntaps, [0, q, 1.08*q, 0.5], [1,0])
return b

def genhpf(ntaps, bottomf, Ts):
q = bottomf/(1/Ts)
b = scipy.signal.remez(ntaps, [0, q, 1.08*q, 0.5], [0,1])
return b

def genbandpassf(ntaps, bottomf, topf, Ts):
q = bottomf/(1/Ts)
r = topf/(1/Ts)

b = scipy.signal.remez(ntaps, [0,q*0.98,q,r,r*1.02,0.5], [0,1,0])
return b

In [23]:
Ts = 1.0/10000.0

plotfilter( genhpf(100,  500.0, Ts), Ts)
plotfilter( genlpf(100, 3000.0, Ts), Ts)
plotfilter( genbandpassf(100, 1500, 2500, Ts), Ts)


3.9. As 3.8, but with Ts=1/20000.0

In [26]:
Ts = 1.0/20000.0

plotfilter( genhpf(100,  500.0, Ts), Ts)
plotfilter( genlpf(100, 3000.0, Ts), Ts)
plotfilter( genbandpassf(100, 1500, 2500, Ts), Ts)


3.10. Let x1 be a cos of f=800, x2(t) of f=2000 and x3 of f=4500. Let x(t) = x1 + 0.5x2 + 2x3

In [27]:
def plotfilter2(b, Ts):
figure()
time = 3.0
tp = 2*pi*linspace(0, time, time/Ts)
x = cos(800*tp) + 0.5*cos(2000*tp) + 2*cos(4500*tp)
y = scipy.signal.lfilter(b,1,x)
plotspec(y, Ts)

plotfilter2( genhpf(100,  500.0, Ts), Ts)
plotfilter2( genlpf(100, 3000.0, Ts), Ts)
plotfilter2( genbandpassf(100, 1500, 2500, Ts), Ts)