#
# 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)
# 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()
%==============================================================================
%
% Gray Code Generator for Verilog
%
% This script generates a Verilog header (.vh) file which contains a set of
% of `define directives to establish Gray-coded states for a state machine.
% The output of this script is intended to be used in conjunction with a
% Verilog design which uses Gray coding in some application. For more
% information on Gray coding, please see
% http://mathworld.wolfram.com/GrayCode.html.
%
% Usage:
%
% The bit width of the Gray coded states is set as a script variable.
% 2^bit_width states are generated. The output filename is also
% set as a script variable.
%
% Author: Tom Chatt
%
%------------------------------------------------------------------------------
% Revision History:
%
% Original Version Tom Chatt 11/2/2010
%
%==============================================================================
clear all; clc;
%------------------------------------------------------------------------------
bit_width = 8; % Bit width of the generated code
% Print the generated code entries to a Verilog defines file?
print_to_v_def_file = true;
output_file = ['gray_code_',num2str(bit_width),'.vh']; % Output file name
% Gray code state name prefix. State names will be of the form
% GC<bit_width>B_<state number>.
prefix = ['GC',num2str(bit_width),'B_'];
%--------------Do not edit code below this point-------------------------------
% Generate the code, starting with 1 bit code and looping upward
code = [0; 1];
if bit_width ~= 1
for i = 2 : bit_width
code = [zeros(size(code,1), 1), code; ones(size(code,1), 1), flipud(code)]; %#ok<AGROW>
end
end
% Print out the generated code
fid = fopen(output_file, 'w');
for i = 1 : size(code,1)
macro_name = [prefix, num2str(i)];
macro_val = [num2str(bit_width),'''', dec2bin(binvec2dec(code(i,:)),bit_width)];
fprintf(fid, ['ifndef ',macro_name,'\n']);
fprintf(fid, [' ','`define ', macro_name, ' ', macro_val, '\n']);
fprintf(fid, 'endif\n');
end
fclose('all');
# 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()
#
# 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)
%==============================================================================
%
% Gray Code Generator for Verilog
%
% This script generates a Verilog header (.vh) file which contains a set of
% of `define directives to establish Gray-coded states for a state machine.
% The output of this script is intended to be used in conjunction with a
% Verilog design which uses Gray coding in some application. For more
% information on Gray coding, please see
% http://mathworld.wolfram.com/GrayCode.html.
%
% Usage:
%
% The bit width of the Gray coded states is set as a script variable.
% 2^bit_width states are generated. The output filename is also
% set as a script variable.
%
% Author: Tom Chatt
%
%------------------------------------------------------------------------------
% Revision History:
%
% Original Version Tom Chatt 11/2/2010
%
%==============================================================================
clear all; clc;
%------------------------------------------------------------------------------
bit_width = 8; % Bit width of the generated code
% Print the generated code entries to a Verilog defines file?
print_to_v_def_file = true;
output_file = ['gray_code_',num2str(bit_width),'.vh']; % Output file name
% Gray code state name prefix. State names will be of the form
% GC<bit_width>B_<state number>.
prefix = ['GC',num2str(bit_width),'B_'];
%--------------Do not edit code below this point-------------------------------
% Generate the code, starting with 1 bit code and looping upward
code = [0; 1];
if bit_width ~= 1
for i = 2 : bit_width
code = [zeros(size(code,1), 1), code; ones(size(code,1), 1), flipud(code)]; %#ok<AGROW>
end
end
% Print out the generated code
fid = fopen(output_file, 'w');
for i = 1 : size(code,1)
macro_name = [prefix, num2str(i)];
macro_val = [num2str(bit_width),'''', dec2bin(binvec2dec(code(i,:)),bit_width)];
fprintf(fid, ['ifndef ',macro_name,'\n']);
fprintf(fid, [' ','`define ', macro_name, ' ', macro_val, '\n']);
fprintf(fid, 'endif\n');
end
fclose('all');