DSPRelated.com
   Code Snippets

Halfband Filter Design with Python/Scipy

Christopher Felton ● February 24, 2012●5 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 = fig.add_subplot(211)
ax0.stem(numpy.arange(len(h)), h)
ax0.grid(True)
ax0.set_title('Parks-McClellan (remez) Impulse Response')
ax1 = fig.add_subplot(212)
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 = fig.add_subplot(111)
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_xlabel('Normalized Frequency (radians)')
ax1.set_title('Half Band Filter Frequency Response')
fig.savefig('hb_rsp.png')

Python zplane function

Christopher Felton ● December 17, 2011●4 comments ● Coded in Python 
#
# Copyright (c) 2011 Christopher Felton
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# 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
    uc = patches.Circle((0,0), radius=1, fill=False,
                        color='black', ls='dashed')
    ax.add_patch(uc)

    # 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

Christopher Felton ● 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()

Fast audio convolution

Hernán Ordiales ● November 19, 2010 ● Coded in Python 
#!/usr/bin/python
# -*- coding: utf-8 -*-

#
# Copyright (C) 2006 by Hernán Ordiales
# <audiocode@uint8.com.ar>
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#

# Fast convolution of an audio signal with the impulse response of a system modeled as LTI (Linear Time Invariant) one.
# Use: fast_conv.py source.wav impulse_response.wav output.wav

from wav_array import *
from FFT import fft, inverse_fft
from pylab import size
import sys

def nextpow2( L ):
	N = 2
	while N < L: N = N * 2
	return N
	
def fast_conv_vect( x, h ):
	# searches for the amount of points required to perform the FFT
	L = size(h) + size(x) - 1 # linear convolution length
	N = nextpow2( L )
	# Note: N>=L is needed because the IDFT of the multiplication is the circular convolution and to match it to the
	# common one, N>=L is required (where L=N1+N2-1;N1=length(x);N2=length(h))

	# FFT(X,N) is the N points FFT, with zero padding if X has less than N points and truncated if has more.
	H = fft( h, N ) # Fourier transform of the impulse
	X = fft( x, N ) # Fourier transform of the input signal
	
	Y = H * X # spectral multiplication
	y = inverse_fft( Y ) # time domain again
	return y

source = sys.argv[1]
impulse = sys.argv[2]
conv_output = sys.argv[3]

clip_factor = 1.01 # clip factor default value

[ h1, Fs1, h_bits ] = wavread( impulse ) # impulse response
[ x1, Fs2, x_bits ] = wavread( source ) # file to process

if Fs1 == Fs2 : # if sample rates are the same
	print "Processing..."
	y1 = fast_conv_vect( x1, h1 ).real # takes the real part to avoid a too small complex part (around e-18)

	# audio normalization: if "y = y/max(y)" -> "values out of [-1,+1] range are clipped"
	y1 = y1/( max(y1)*clip_factor ) # to avoid clipping
	
	wavwrite( y1, Fs1, conv_output )

	print "Output file:", conv_output
	print "Convolution success."
else:
	print "Error: files has different sample rate."

Recursive FFT in Python Convertible to Verilog/VHDL

Christopher Felton ● November 3, 2010●2 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)