## Halfband Filter Design with Python/Scipy

February 24, 20125 comments Coded in Python
``````# The following is a Python/scipy snippet to generate the
# coefficients for a halfband filter.  A halfband filter
# is a filter where the cutoff frequency is Fs/4 and every
# other coeffecient is zero except the cetner tap.
# Note: every other (even except 0) is 0, most of the coefficients
#       will be close to zero, force to zero actual

import numpy
from numpy import log10, abs, pi
import scipy
from scipy import signal
import matplotlib
import matplotlib.pyplot
import matplotlib as mpl

# ~~[Filter Design with Parks-McClellan Remez]~~
N = 32  # Filter order
# Filter symetric around 0.25 (where .5 is pi or Fs/2)
bands = numpy.array([0., .22, .28, .5])
h = signal.remez(N+1, bands, [1,0], [1,1])
h[abs(h) <= 1e-4] = 0.
(w,H) = signal.freqz(h)

# ~~[Filter Design with Windowed freq]~~
b = signal.firwin(N+1, 0.5)
b[abs(h) <= 1e-4] = 0.
(wb, Hb) = signal.freqz(b)

# Dump the coefficients for comparison and verification
print('          remez       firwin')
print('------------------------------------')
for ii in range(N+1):
print(' tap %2d   %-3.6f    %-3.6f' % (ii, h[ii], b[ii]))

## ~~[Plotting]~~
# Note: the pylab functions can be used to create plots,
#       and these might be easier for beginners or more familiar
#       for Matlab users.  pylab is a wrapper around lower-level
#       MPL artist (pyplot) functions.
fig = mpl.pyplot.figure()
ax0.stem(numpy.arange(len(h)), h)
ax0.grid(True)
ax0.set_title('Parks-McClellan (remez) Impulse Response')
ax1.stem(numpy.arange(len(b)), b)
ax1.set_title('Windowed Frequency Sampling (firwin) Impulse Response')
ax1.grid(True)
fig.savefig('hb_imp.png')

fig = mpl.pyplot.figure()
ax1.plot(w, 20*log10(abs(H)))
ax1.plot(w, 20*log10(abs(Hb)))
ax1.legend(['remez', 'firwin'])
bx = bands*2*pi
ax1.axvspan(bx[1], bx[2], facecolor='0.5', alpha='0.33')
ax1.plot(pi/2, -6, 'go')
ax1.axvline(pi/2, color='g', linestyle='--')
ax1.axis([0,pi,-64,3])
ax1.grid('on')
ax1.set_ylabel('Magnitude (dB)')
ax1.set_title('Half Band Filter Frequency Response')
fig.savefig('hb_rsp.png')``````

## Python zplane function

December 17, 20114 comments Coded in Python
``````#
# Copyright (c) 2011 Christopher Felton
#
# This program is free software: you can redistribute it and/or modify
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
#

# The following is derived from the slides presented by
# Alexander Kain for CS506/606 "Special Topics: Speech Signal Processing"
# CSLU / OHSU, Spring Term 2011.

import numpy as np
import matplotlib.pyplot as plt
from  matplotlib import patches
from matplotlib.figure import Figure
from matplotlib import rcParams

def zplane(b,a,filename=None):
"""Plot the complex z-plane given a transfer function.
"""

# get a figure/plot
ax = plt.subplot(111)

# create the unit circle
color='black', ls='dashed')

# The coefficients are less than 1, normalize the coeficients
if np.max(b) > 1:
kn = np.max(b)
b = b/float(kn)
else:
kn = 1

if np.max(a) > 1:
kd = np.max(a)
a = a/float(kd)
else:
kd = 1

# Get the poles and zeros
p = np.roots(a)
z = np.roots(b)
k = kn/float(kd)

# Plot the zeros and set marker properties
t1 = plt.plot(z.real, z.imag, 'go', ms=10)
plt.setp( t1, markersize=10.0, markeredgewidth=1.0,
markeredgecolor='k', markerfacecolor='g')

# Plot the poles and set marker properties
t2 = plt.plot(p.real, p.imag, 'rx', ms=10)
plt.setp( t2, markersize=12.0, markeredgewidth=3.0,
markeredgecolor='r', markerfacecolor='r')

ax.spines['left'].set_position('center')
ax.spines['bottom'].set_position('center')
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)

# set the ticks
r = 1.5; plt.axis('scaled'); plt.axis([-r, r, -r, r])
ticks = [-1, -.5, .5, 1]; plt.xticks(ticks); plt.yticks(ticks)

if filename is None:
plt.show()
else:
plt.savefig(filename)

return z, p, k``````

## FPGA IIR Lowpass Direct Form I Filter Generator

August 24, 2011 Coded in Python
``````# Instantiate the SIIR object.  Pass the cutoff frequency
# Fc and the sample rate Fs in Hz.  Also define the input
# and output fixed-point type.  W=(wl, iwl) where
# wl = word-length and iwl = integer word-length.  This
# example uses 23 fraction bits and 1 sign bit.
>>> from siir import SIIR
>>> flt = SIIR(Fstop=1333, Fs=48000, W=(24,0))

# Plot the response of the fixed-point coefficients
>>> plot(flt.hz, 20*log10(flt.h)

# Create a testbench and run a simulation
# (get the simulated response)
>>> from myhdl import Simulation
>>> tb = flt.TestFreqResponse(Nloops=128, Nfft=1024)
>>> Simulation(tb).run()
>>> flt.PlotResponse()

# Happy with results generate the Verilog and VHDL
>>> flt.Convert()``````

## Recursive FFT in Python Convertible to Verilog/VHDL

November 3, 20102 comments Coded in Python
``````#
# rfft.py
#
# This file contains a recursive version of the fast-fourier transform and
# support test functions.  This module utilizes the numpy (numpy.scipy.org)
# library.
#
# References
#   - http://www.cse.uiuc.edu/iem/fft/rcrsvfft/
#   - "A Simple and Efficient FFT Implementation in C++", by Vlodymyr Myrnyy

import numpy
from numpy.fft import fft
from numpy import sin, cos, pi, ones, zeros, arange, r_, sqrt, mean

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
def rFFT(x):
"""
Recursive FFT implementation.

References
-- http://www.cse.uiuc.edu/iem/fft/rcrsvfft/
-- "A Simple and Efficient FFT Implementation in C++"
by Vlodymyr Myrnyy
"""

n = len(x)

if (n == 1):
return x

w = getTwiddle(n)
m = n/2;
X = ones(m, float)*1j
Y = ones(m, float)*1j

for k in range(m):
X[k] = x[2*k]
Y[k] = x[2*k + 1]

X = rFFT(X)
Y = rFFT(Y)

F = ones(n, float)*1j
for k in range(n):
i = (k%m)
F[k] = X[i] + w[k] * Y[i]

return F

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
def getTwiddle(NFFT=8):
"""Generate the twiddle factors"""

W = r_[[1.0 + 1.0j]*NFFT]

for k in range(NFFT):
W[k] = cos(2.0*pi*k/NFFT) - 1.0j*sin(2.0*pi*k/NFFT)

return W

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
def DFT(x, N=8):
"""
Use the direct definition of DFT for verification
"""
y = [1.0 + 1.0j]*N
y = r_[y]
for n in range(N):
wsum = 0 + 0j;
for k in range(N):
wsum = wsum + (cos(2*pi*k*n/N) - (1.0j * sin(2*pi*k*n/N)))*x[k]

y[n] = wsum

return y

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
def test_rfft(N      = 64,      # FFT order to test
nStart = 0.2,     # Note aliased signal included
nStep  = 2.1,     # Samples per period step
pStep  = pi/4,    # Phase step size
limErr = 10e-12,  # Error limit to check
maxErr = 0        # Max difference
):
"""
Use the built in numpy FFT functions and the direct
implemenation of the DFT to verify the recursive FFT.

This testbench verifies the different implementations are within
a certain limit.  Because of the different implemenations the values
could be slightly off (computer representation calculation error).
"""

# Use test signal nStart:nStep:N samples per cycle
for s in arange(nStart, N+nStep, nStep):
for p in arange(0, pi+pStep, pStep):

n = arange(N, 0, -1)
x = cos(2*pi*n/s + p)

xDFT = DFT(x,N)
nFFT = fft(x,N)
xFFT = rFFT(x)

rmsErrD = sqrt(mean(abs(xDFT - xFFT))**2)
rmsErrN = sqrt(mean(abs(nFFT - xFFT))**2)

if rmsErrD > limErr or rmsErrN > limErr:
print s, p, "Error!", rmsErrD, rmsErrN
print xDFT
print nFFT
print xFFT

if rmsErrD > maxErr:
maxErr = rmsErrD
elif rmsErrN > maxErr:
maxErr = rmsErrN

print "N %d maxErr = %f " % (N,maxErr)

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# If the module is run test a bunch of different size FFTs
if __name__ == '__main__':

# The following is fairly exhaustive and will take some time
# to run.
tv = 2**arange(1,12)
for nfft in tv:
test_rfft(N=nfft)``````