Commit dc2c0e65 by Javier D. Garcia-Lasheras

### Cleaning the tree

parent 963d8c15
This diff is collapsed.
 # This file is part of librefdatool. librefdatool 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, version 2. # # 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 General Public License for more # details. # # You should have received a copy of the GNU General Public License along with # this program; if not, write to the Free Software Foundation, Inc., 51 # Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. # # Copyright (C) 2013 Javier D. Garcia-Lasheras # # Analysis class code is based in previous open source work: # * mfreqz: by Matti Pastell # * zplane: by Chris Felton import numpy as np from numpy import pi, log10 from scipy import signal from matplotlib import pyplot as plt from matplotlib import mlab from matplotlib import patches from matplotlib.figure import Figure from matplotlib import rcParams class Analysis(): '''This class include tools for LTI filter analysis, sporting math models for studying coefficient quantization effects. Note (limited to FIR filters) ''' def analyze_zero_pole(self, c, p, q, filename=None): '''Plot graphical zero/pole analysis in Z-plane: * c: FIR filter coefficients array. * filename: optional file for storing plot and not showing it. ''' print('Plotting Z Analysis...') # Quantize coefficients d = np.zeros(len(c)) for ii in range(len(c)): d[ii] = int(c[ii]*(2**(p+q-1))) # Temporal assignation: only valid for FIR print('Coefficients') b1 = c a1 = 1 print(b1) print(np.max(b1)) b2 = d a2 = 1 print(b2) print(np.max(b1)) plt.figure('libre-fdatool analysis'); # 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) # 1 - The coefficients are less than 1, normalize the coeficients if np.max(b1) > 1: kn1 = np.max(b1) b1 = b1/float(kn1) else: kn1 = 1 if np.max(a1) > 1: kd1 = np.max(a1) a1 = a1/float(kd1) else: kd1 = 1 # 2 - The coefficients are less than 1, normalize the coeficients if np.max(b2) > 1: kn2 = np.max(b2) b2 = b2/float(kn2) else: kn2 = 1 if np.max(a2) > 1: kd2 = np.max(a2) a2 = a2/float(kd2) else: kd2 = 1 # Get the poles and zeros p1 = np.roots(a1) z1 = np.roots(b1) k1 = kn1/float(kd1) p2 = np.roots(a2) z2 = np.roots(b2) k2 = kn2/float(kd2) # 1 - Plot the zeros and set marker properties tz1 = plt.plot(z1.real, z1.imag, 'bo', ms=10) plt.setp( tz1, markersize=10.0, markeredgewidth=1.0, markeredgecolor='b', markerfacecolor='b') # 1 - Plot the poles and set marker properties tp1 = plt.plot(p1.real, p1.imag, 'bx', ms=10) plt.setp( tp1, markersize=12.0, markeredgewidth=3.0, markeredgecolor='b', markerfacecolor='b') # 2 - Plot the zeros and set marker properties tz2 = plt.plot(z2.real, z2.imag, 'ro', ms=10) plt.setp( tz2, markersize=10.0, markeredgewidth=1.0, markeredgecolor='r', markerfacecolor='r') # 2 - Plot the poles and set marker properties tp2 = plt.plot(p2.real, p2.imag, 'rx', ms=10) plt.setp( tp2, markersize=12.0, markeredgewidth=3.0, markeredgecolor='r', markerfacecolor='r') # set axis plt.axis('scaled') plt.ylabel('Imaginary Component') plt.xlabel('Real Component') plt.title(r'Zero-Pole Diagram (Blue=Float; Red=Int)') plt.grid(True) # show or store plot return plt def analyze_frequency_response(self, c, p, q, filename=None): '''Plot graphical Magnitude/Phase analysis in frequency domain: * c: FIR filter coefficients array. * filename: optional file for storing plot and not showing it. ''' # Quantize coefficients d = np.zeros(len(c)) for ii in range(len(c)): d[ii] = int(c[ii]*(2**(p+q-1))) d = d/(2**(p+q-1)) print(c) print(d) wc,hc = signal.freqz(c,1) hc_dB = 20 * log10 (abs(hc)) wd,hd = signal.freqz(d,1) hd_dB = 20 * log10 (abs(hd)) plt.figure('libre-fdatool analysis'); plt.subplot(211) plt.plot(wc/max(wc),hc_dB,'b') plt.plot(wd/max(wd),hd_dB,'r') plt.ylim(-150, 5) plt.ylabel('Magnitude (dB)') plt.xlabel(r'Normalized Frequency (x$\pi$rad/sample)') plt.title(r'Magnitude response (Blue=Float; Red=Int)') plt.grid(True) plt.subplot(212) hc_Phase = np.unwrap(np.arctan2(np.imag(hc),np.real(hc))) hd_Phase = np.unwrap(np.arctan2(np.imag(hd),np.real(hd))) plt.plot(wc/max(wc),hc_Phase,'b') plt.plot(wd/max(wd),hd_Phase,'r') plt.ylabel('Phase (radians)') plt.xlabel(r'Normalized Frequency (x$\pi$rad/sample)') plt.title(r'Phase response (Blue=Float; Red=Int)') plt.grid(True) plt.subplots_adjust(hspace=0.5) return plt
This source diff could not be displayed because it is too large. You can view the blob instead.
 # This file is part of librefdatool. librefdatool 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, version 2. # # 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 General Public License for more # details. # # You should have received a copy of the GNU General Public License along with # this program; if not, write to the Free Software Foundation, Inc., 51 # Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. # # Copyright (C) 2013 Javier D. Garcia-Lasheras # import numpy as np from numpy import pi, log10 from scipy import signal from matplotlib import pyplot as plt from matplotlib import mlab class control: '''This class acts as control interface for filter parameters, such as structure realization, fixed-point approximation... Note (limited to FIR filters) ''' def __init__(self): '''This method provides initial filter parameters, executing each set command with default values when a new filter object is intantiated. ''' self.setBus() self.setScaling() self.setStructure() self.setModel() self.setTrace() self.setWorkspace() self.setToolchain() def setBus(self, busX=[1,15], busY=[1,15], busC=[1,15]): '''assign the different bus width: The value is given as a pair of natural numbers [p,q]: - p is the number of bits representing the integer part. - q is the number of bits representing the fractional part There are three different configurable buses: * busX: Input bit width * busY: Output bit width * busC: Coefficient bit width ''' self.busX = busX self.busY = busY self.busC = busC def setScaling(self, scalingX=16, scalingC=16): '''assign the different FP scaling factor: * scalingX: Input FP scaling factor (power of two exponent) * scalingC: Coefficient FP scaling factor (power of two exponent) ''' self.scalingX = scalingX self.scalingC = scalingC def setStructure(self, structure='firfilt'): '''assign filter structure or realization. Availability depends on python-to-HDL toolchain selection. TBD: only "firfilt" structure is available Future structures are type I/II direct forms, transposed... ''' self.structure = structure def setTrace(self, trace=False): '''assign activation state for VCD file tracer. If true, a *.vcd file will be generated by simulation process. This file can be analyzed with third party tools such as GTKWave. ''' self.trace = trace def setModel(self, model='myhdl'): '''assign HDL model, this is, the languaje for the toolchain. When not selected the native/default value, it relies on languaje conversion cogeneration & Icarus cosimulation. For myhdl toolchain, valid values are [myhdl, verilog, vhdl, vhdl2] TBD: this mechanism still has to be adapted for toolchain selection and including migen an other tools (maybe GHDL for cosimulation). ''' self.model = model def setWorkspace(self, workspace='.'): '''assign workspace folder for filter design and analysis. By default, the path is the one where the librefdatool is called. In this path, output files will be generated: * VHDL, Verilog converted files * VCD simulation files ''' self.workspace = workspace def setToolchain(self, toolchain='myhdl'): '''Choose python-to-HDL toolchain. By default, MyHDL is the toolchain in use. TBD: migen toolchain will be developed in the near future ''' self.toolchain = toolchain
core.py deleted 100644 → 0
 # This file is part of librefdatool. librefdatool 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, version 2. # # 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 General Public License for more # details. # # You should have received a copy of the GNU General Public License along with # this program; if not, write to the Free Software Foundation, Inc., 51 # Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. # # Copyright (C) 2013 Javier D. Garcia-Lasheras # import sys import os from simcore import simcore class core: '''This class provides the principal librefdatool methods, including running signal filtering with HDL models. TBD: limited to FIR filter simulation ''' def librefilter(self, b, a, x): '''This method provides filter hardware simulation, using a HDL description & a simulation engine parametrized with the device control values. * c: FIR filter coefficients array. * x: Input signal array. * return: simulated output signal NOTE: the filter behaviour is tuned to match with: scipy.signal.lfilter (the purpose is direct comparision of outputs) ''' # Save path and jump to workspace (TBD: if exists!!) self.savedPath = os.getcwd() os.chdir(self.workspace) runner = simcore(); response = runner.simfilter(b, a, x, self.structure, self.model, self.trace, self.busX, self.busY, self.busC, self.scalingX, self.scalingC) # Return from workspace folder os.chdir(self.savedPath) return response
device.py deleted 100644 → 0
 # This file is part of librefdatool. librefdatool 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, version 2. # # 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 General Public License for more # details. # # You should have received a copy of the GNU General Public License along with # this program; if not, write to the Free Software Foundation, Inc., 51 # Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. # # Copyright (C) 2013 Javier D. Garcia-Lasheras # # LibreFDATool Libraries from core import * from analysis import * from control import * class device( control, core, Analysis): '''This class instantiate a librefdatool filter device object. A device represents an HDL model that can be widely parametrized and used in python signal processing analysis and simulation. The device class act as a main wrapper that includes function subclasses TBD: only FIR filter is supported. General LTI system is in development. ''' pass
filter.py deleted 100644 → 0
This diff is collapsed.
This diff is collapsed.

2.93 KB

3.12 KB

3.05 KB

This diff is collapsed.
scope.py deleted 100644 → 0
 # This file is part of librefdatool. librefdatool 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, version 2. # # 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 General Public License for more # details. # # You should have received a copy of the GNU General Public License along with # this program; if not, write to the Free Software Foundation, Inc., 51 # Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. # # Copyright (C) 2013 Javier D. Garcia-Lasheras # import warnings import numpy as np from numpy import pi, log10 from scipy import signal from matplotlib import pyplot as plt from matplotlib import patches, mlab import sys, os, random '''This module include the figure plotters used in Libre-FDATool. ''' class Waveform: def __init__(self, value=[0], label='waveform', color='#000000', linestyle='-', marker='o', markersize=8): self.value = value self.label = label self.color = color self.linestyle = linestyle self.marker = marker self.markersize = markersize def analyze_pole_zero(figure, b, a, p, q, grid): '''Plot graphical zero/pole analysis in Z-plane: * b: LTI transfer function Numerator. * a: LTI transfer function Denominator. * filename: optional file for storing plot and not showing it. ''' print('Plotting Z Analysis...') #TODO: use DAC / ADC like function, this is bad!!! # Quantize coefficients b2 = np.zeros(len(b)) for ii in range(len(b)): b2[ii] = int(b[ii]*(2**(p+q-1))) b2 = b2/(2**(p+q-1)) a2 = np.zeros(len(a)) for ii in range(len(a)): a2[ii] = int(a[ii]*(2**(p+q-1))) a2 = a2/(2**(p+q-1)) # Temporal assignation: only valid for FIR print('Coefficients') b1 = b a1 = a # 1 - The coefficients must be less than 1, normalize the coefficients if np.max(b1) > 1: kn1 = np.max(b1) b1 = b1/float(kn1) else: kn1 = 1 if np.max(a1) > 1: kd1 = np.max(a1) a1 = a1/float(kd1) else: kd1 = 1 # 2 - The coefficients must be than 1, normalize the coefficients if np.max(b2) > 1: kn2 = np.max(b2) b2 = b2/float(kn2) else: kn2 = 1 if np.max(a2) > 1: kd2 = np.max(a2) a2 = a2/float(kd2) else: kd2 = 1 warningMessage = '' warnings.simplefilter('error') try: # Get the poles and zeros p1 = np.roots(a1) z1 = np.roots(b1) p2 = np.roots(a2) z2 = np.roots(b2) except ValueError as exVE: warningMessage = '%s' % exVE except ZeroDivisionError as exZDE: warningMessage = '%s' % exZDE except RuntimeWarning as exRW: warningMessage = '%s' % exRW if warningMessage != '': return False, warningMessage else: # Clear the figure figure.clear() # get a figure/plot poleZeroSubplot = figure.add_subplot(111) # create the unit circle uc = patches.Circle((0,0), radius=1, fill=False, color='black', ls='dashed') #ax.add_patch(uc) poleZeroSubplot.add_patch(uc) # 1 - Plot the zeros and set marker properties poleZeroSubplot.plot(z1.real, z1.imag, 'bo', ms=10) # 1 - Plot the poles and set marker properties poleZeroSubplot.plot(p1.real, p1.imag, 'bx', ms=10) # 2 - Plot the zeros and set marker properties poleZeroSubplot.plot(z2.real, z2.imag, 'ro', ms=10) # 2 - Plot the poles and set marker properties poleZeroSubplot.plot(p2.real, p2.imag, 'rx', ms=10) # set axis poleZeroSubplot.axis('scaled') poleZeroSubplot.set_ylabel('Imaginary Component') poleZeroSubplot.set_xlabel('Real Component') poleZeroSubplot.set_title(r'Zero-Pole Diagram (Blue=Float; Red=Int)') poleZeroSubplot.grid(grid) return True, 'OK' def analyze_frequency_response(figure, b, a, p, q, grid): '''Plot graphical Magnitude/Phase analysis in frequency domain: * c: FIR filter coefficients array. * filename: optional file for storing plot and not showing it. ''' # Quantize coefficients b2 = np.zeros(len(b)) for ii in range(len(b)): b2[ii] = int(b[ii]*(2**(p+q-1))) b2 = b2/(2**(p+q-1)) a2 = np.zeros(len(a)) for ii in range(len(a)): a2[ii] = int(a[ii]*(2**(p+q-1))) a2 = a2/(2**(p+q-1)) # TODO: If the coefficients are extremely low, the discretized version # may be equal to zero, which suppose a log10 crash - divide by zero! # We need to rise an advice about rising the coefficient bits!! warningMessage = '' warnings.simplefilter('error') try: wc,hc = signal.freqz(b,a) hc_dB = 20 * log10 (abs(hc)) wd,hd = signal.freqz(b2,a2) hd_dB = 20 * log10 (abs(hd)) except ValueError as exVE: warningMessage = '%s' % exVE except ZeroDivisionError as exZDE: warningMessage = '%s' % exZDE except RuntimeWarning as exRW: warningMessage = '%s' % exRW figure.clear() magnitudeSubplot = figure.add_subplot(211) magnitudeSubplot.set_ylim(-150, 5) magnitudeSubplot.set_ylabel('Magnitude (dB)') magnitudeSubplot.set_xlabel(r'Normalized Frequency (x$\pi$rad/sample)') magnitudeSubplot.grid(grid) phaseSubplot = figure.add_subplot(212) phaseSubplot.set_ylabel('Phase (radians)') phaseSubplot.set_xlabel(r'Normalized Frequency (x$\pi$rad/sample)') phaseSubplot.grid(grid) figure.subplots_adjust(hspace=0.5) if warningMessage != '': magnitudeSubplot.plot(wc/max(wc),hc_dB,'b') hc_Phase = np.unwrap(np.arctan2(np.imag(hc),np.real(hc))) magnitudeSubplot.set_title(r'Magnitude response (Float=Blue; Int=Error)') phaseSubplot.plot(wc/max(wc),hc_Phase,'b') phaseSubplot.set_title(r'Phase response (Float=Blue; Int=Error)') return False, warningMessage else: magnitudeSubplot.plot(wc/max(wc),hc_dB,'b') magnitudeSubplot.plot(wd/max(wd),hd_dB,'r') magnitudeSubplot.set_title(r'Magnitude response (Float=Blue; Int=Red)') hc_Phase = np.unwrap(np.arctan2(np.imag(hc),np.real(hc))) hd_Phase = np.unwrap(np.arctan2(np.imag(hd),np.real(hd))) phaseSubplot.plot(wc/max(wc),hc_Phase,'b') phaseSubplot.plot(wd/max(wd),hd_Phase,'r') phaseSubplot.set_title(r'Phase response (Float=Blue; Int=Error)') return True, 'OK' def scopeTime(figure, waveform, grid): '''this method shows a dual time domain scope: * s1: channel-1 input signal (blue, float) * s2: channel-2 input signal (red, int) ''' # Simulate float system response: print('Run Scope Time') figure.clear() scopeTimeSubplot = figure.add_subplot(111) scopeTimeSubplot.set_ylabel('Value') scopeTimeSubplot.set_xlabel('Sample') scopeTimeSubplot.set_title(r'Time Scope') scopeTimeSubplot.grid(grid) for ii in range(len(waveform)): scopeTimeSubplot.plot(waveform[ii].value, color = waveform[ii].color, label = waveform[ii].label, linestyle = waveform[ii].linestyle, marker = waveform[ii].marker, markersize = waveform[ii].markersize) scopeTimeSubplot.legend().draggable(state=True, use_blit=True) def scopePower(figure, waveform, grid): '''this method shows the estimated power spectrum for a signal: * s: channel-1 assigned signal (blue) ''' print('Plotting Power Spectrum...') figure.clear() scopePowerSubplot = figure.add_subplot(111) scopePowerSubplot.set_ylabel('Power (dB)') scopePowerSubplot.set_xlabel('Normalized Frequency') scopePowerSubplot.set_title(r'Signal Power') scopePowerSubplot.grid(grid) for ii in range(len(waveform)): Ps,fs = mlab.psd(waveform[ii].value) scopePowerSubplot.plot(fs, 10*log10(abs(Ps)), color = waveform[ii].color, label = waveform[ii].label, linestyle = waveform[ii].linestyle) scopePowerSubplot.legend().draggable(state=True, use_blit=True) def scopeError(figure, s1, s2, grid): '''this method shows an error analysis between input signals: * s1: channel-1 input signal * s1: channel-2 input signal ''' # *** Error Analisys *** sdiff = np.abs(s1 - s2) print('- Maximum error = ', np.max(sdiff)) print('- Mean error = ', np.mean(sdiff**2)) # Check for error tolerance # assert np.max(sdiff) < 1e-3, "check if error is too large" # Plot Error report figure.clear() scopeErrorSublot = figure.add_subplot(111)