Plots the single-bin frequency response of direct FFT and PFB.

singlebinresponse.py

#!/usr/bin/python

# singlebinresponse.py
# Plots the single-bin frequency response of direct FFT and PFB
#
# Created by Jayanth Chennamangalam
# (PFB coefficient generation based on code by Sean McHugh, UCSB)

import sys
import getopt
import numpy as np
import matplotlib.pyplot as plt


# function definitions
def PrintUsage(ProgName):
    "Prints usage information."
    print("Usage: " + ProgName + " [options]")
    print("    -h  --help                 Display this usage information")
    print("    -n  --nfft <value>         Number of points in FFT")
    print("    -t  --taps <value>         Number of taps in PFB")
    return


def safeLog10(x, minval=0.0000000001):
    return np.log10(x.clip(min=minval))


# default values
NFFT = 1024                 # number of points in FFT
NTaps = 8                   # number of taps in PFB

# get the command line arguments
ProgName = sys.argv[0]
OptsShort = "hn:t:"
OptsLong = ["help", "nfft=", "taps="]

# get the arguments using the getopt module
try:
    (Opts, Args) = getopt.getopt(sys.argv[1:], OptsShort, OptsLong)
except getopt.GetoptError as ErrMsg:
    # print usage information and exit
    sys.stderr.write("ERROR: " + str(ErrMsg) + "!\n")
    PrintUsage(ProgName)
    sys.exit(1)

# parse the arguments
for o, a in Opts:
    if o in ("-h", "--help"):
        PrintUsage(ProgName)
        sys.exit()
    elif o in ("-n", "--nfft"):
        NFFT = int(a)
    elif o in ("-t", "--taps"):
        NTaps = int(a)
    else:
        PrintUsage(ProgName)
        sys.exit(1)

M = NTaps * NFFT

# generate PFB coefficients
X = np.array([(float(i) / NFFT) - (float(NTaps) / 2) for i in range(M)])
PFBCoeff = np.sinc(X) * np.hanning(M)

# take FFT of the PFB coefficients - make sure N is large enough compared to
# NFFT * NTaps
N = 1048576
f = np.arange(-float(NFFT / 2), float(NFFT / 2), float(NFFT) / N)
PFBSpec = np.fft.fft(PFBCoeff, n=N)

# generate 'no-PFB' coefficients, i.e., a rectangular window (direct FFT)
NoPFBCoeff = np.ones(NFFT)
# take FFT of the rectangular window
NoPFBSpec = np.fft.fft(NoPFBCoeff, n=N)

# plot the frequency responses of the PFB and direct FFT
SpecPlot = 20 * safeLog10(np.abs(np.fft.fftshift(PFBSpec)))
NoPFBSpecPlot = 20 * safeLog10(np.abs(np.fft.fftshift(NoPFBSpec)))
SpecPlot = SpecPlot - np.max(SpecPlot)
NoPFBSpecPlot = NoPFBSpecPlot - np.max(NoPFBSpecPlot)
plt.plot(f, SpecPlot, 'r-')
plt.plot(f, NoPFBSpecPlot, 'b--')
plt.xlabel("Frequency (bins)")
plt.ylabel("Power (dB)")
plt.xlim(-4, 4)
plt.ylim(-100, 0)
plt.legend(('PFB', 'FFT'), 'upper right')
plt.show()