DSP HELP! Python (or C++) options...

Started by will20fitz, October 20, 2017, 06:20:05 AM

Previous topic - Next topic

will20fitz

Hey there, i'm dipping my toes into the DSP water after a number of analog builds.

I have some code in Python (and the potential to run in C++), and was wondering what my options were with regard to DSP chips.

FV-1 has its own language as do some other chips i've researched... Is there something out there that i could use?

Basically i have a Python program that i want to use to output float (or string/tuple) numbers that would form variables within a DSP program.

Is this possible at all?! :icon_eek:

Digital Larry

I'm not sure that Python, being interpreted, is going to give you required performance for doing real time DSP.

You should check out:  https://ccrma.stanford.edu/software/stk/

I've built up a few things this way on Mac OS.  As far as portability to smaller systems, I cannot tell you.
Digital Larry
Want to quickly design your own effects patches for the Spin FV-1 DSP chip?
https://github.com/HolyCityAudio/SpinCAD-Designer

orbitbot

There's a free course on Coursera right now that uses Python and this toolset for audio processing: https://github.com/MTG/sms-tools -- I've been falling behind a bit in the course, but I don't think that much of it is live audio processing.

MicroPython is a variant that runs on the ESP8266 chipset, but I haven't looked into it much. At any rate, there's a newer version of those boards with better perf, but those boards do not have audio quality ADC and DACs as far as I'm aware, so by default you might end up with something of the quality of arduino DSP pedals ( f.e. PedalShield UNO and the Pedal Pi here https://www.electrosmash.com/ ).

markseel

will20fitz it sounds like you want to use Python as a design tool rather than a real-time DSP platform true?  I'll follow up with posts on how to do that.  But *if* you want to use Python with DSP to test out some ideas you can do it in non-real time by processing WAV files and post-processing / analyzing the output.  If you want to do some real time Python DSP, other than very simple algorithms in pure Python try using Cython.  http://bastibe.de/2012-11-02-real-time-signal-processing-in-python.html

markseel

Here's a Python script to compute coefficients for FIR, IIR and single biquad filters for lowpass, highpass, peaking, shevling etc filter responses.  It outputs coefficients in floating point so if you want to use fixed point you can initialize C/C++ coefficient arrays using literal assignment with floating point to integer conversion using the preprocessor:


// FQ converts floating-point to Q31 fixed-point
// QF converts Q31 fixed-point to floating point

#define FQ(hh) (((hh)<0.0)?((int)((double)(1<<31)*(hh)-0.5)):((int)(((double)(1<<31)-1)*(hh)+0.5)))
#define QF(xx) (((int)(xx)<0)?((double)(int)(xx))/(1<<31):((double)(xx))/((1<<31)-1))

// Biquad coefficients (B0,B1,B2,A1,A2) in Q31 fixed point format
// int my_biquad_coeffs[5] = { FQ(+1.0), FQ(0.0), FQ(0.0), FQ(0.0), FQ(0.0) };


A usage example.  It outputs the filter coefficients and plots the magnitude/phase responses.  Note the the frequency is specified relative to the nyquist frequency.  So if the sampling frequency is 48 kHz and your filter frequency is 960 Hz then Fc = 960/48000 = 0.2.  Gain is specified in dB.

bash$ python filters.py biquad type=lowshelf freq=0.2 Q=1.0 gain=12.0 plot=yes
BiQuad Type="lowshelf" Freq=0.200000 Q=1.000 Gain=12.0 ---> b0,...b2,a1,...a2
+1.759364836060580,+0.053441984675100,+0.302146260339349,+0.738408458480827,-0.269660653244002


Here's the Python script.  I haven't verified that the computations are correct so use with caution :-)


import sys, struct, numpy

import numpy as np
import scipy.signal as dsp
import matplotlib.pyplot as plot

np.seterr( all='ignore' )

def coefficients_string( q, b, a ):

t = ""
if q == "":
for i in range(0,len(b)):
if b[i]/a[0] >= 0: t += "+%3.15f" % (b[i]/a[0])
if b[i]/a[0]  < 0: t += "%3.15f"  % (b[i]/a[0])
if i < len(b)-1 or len(a) > 1: t += ","
for i in range(1,len(a)):
if -a[i]/a[0] >= 0: t += "+%3.15f" % (-a[i]/a[0])
if -a[i]/a[0]  < 0: t += "%3.15f"  % (-a[i]/a[0])
if i < len(a)-1: t += ","
else:
#t += "G=%s(%3.12f) " % gain
for i in range(0,len(b)):
if b[i]/a[0] >= 0: t += "%s(+%3.15f)" % (q,b[i]/a[0])
if b[i]/a[0]  < 0: t += "%s(%3.15f)"  % (q,b[i]/a[0])
if i < len(b)-1 or len(a) > 1: t += ","
for i in range(1,len(a)):
if -a[i]/a[0] >= 0: t += "%s(+%3.15f)" % (q,-a[i]/a[0])
if -a[i]/a[0]  < 0: t += "%s(%3.15f)"  % (q,-a[i]/a[0])
if i < len(a)-1: t += ","

return t

def plot_response( bb, aa, xmin=None, xmax=None, ymin=-120.0, ymax=3.0 ):

    fs = 0.5
    w,h = dsp.freqz(bb,aa)
    h_dB = 20 * np.log10 (abs(h))
    if xmin == None: xmin = min(fs*w/np.pi)
    if xmax == None: xmax = max(fs*w/np.pi)
    if ymin == None: ymin = min(h_dB)
    if ymax == None: ymax = max(h_dB)
    plot.subplot(211)
    plot.plot(fs*w/max(w),h_dB)
    plot.xlim(xmin,xmax)
    plot.ylim(ymin,ymax)
    plot.ylabel('Magnitude (db)')
    plot.xlabel('Frequency (Fs = 1.0)')
    plot.title('Frequency response')
    plot.subplot(212)
    h_Phase = np.unwrap(np.arctan2(np.imag(h),np.real(h))) / np.pi * 180.0
    plot.plot(w/max(w),h_Phase)
    plot.xlim(xmin,xmax)
    plot.ylabel('Phase (degrees)')
    plot.xlabel('Frequency (Fs = 1.0)')
    plot.title('Phase response')
    plot.subplots_adjust(hspace=0.5)
    plot.show()

def _make_biquad_notch( filter_freq, q_factor ):

w0 = 2.0 * np.pi * filter_freq
alpha = np.sin(w0)/(2.0 * q_factor)

b0 = +1.0
b1 = -2.0 * cos(w0)
b2 = +1.0
a0 = +1.0 + alpha
a1 = -2.0 * cos(w0)
a2 = +1.0 - alpha

return (b0,b1,b2,a0,a1,a2)

def _make_biquad_lowpass( filter_freq, q_factor ):

w0 = 2.0 * np.pi * filter_freq
alpha = np.sin(w0)/(2 * q_factor)

b0 = (+1.0 - np.cos(w0)) / 2.0
b1 =  +1.0 - np.cos(w0)
b2 = (+1.0 - np.cos(w0)) / 2.0
a0 = +1.0 + alpha
a1 = -2.0 * np.cos(w0)
a2 = +1.0 - alpha

return (b0,b1,b2,a0,a1,a2)

def _make_biquad_highpass( filter_freq, q_factor ):

w0 = 2.0 * np.pi * filter_freq
alpha = np.sin(w0)/(2 * q_factor)

b0 =  (1.0 + np.cos(w0)) / 2.0
b1 = -(1.0 + np.cos(w0))
b2 =  (1.0 + np.cos(w0)) / 2.0
a0 =  +1.0 + alpha
a1 =  -2.0 * np.cos(w0)
a2 =  +1.0 - alpha

return (b0,b1,b2,a0,a1,a2)

def _make_biquad_allpass( filter_freq, q_factor ):

w0 = 2.0 * np.pi * filter_freq
alpha = np.sin(w0)/(2.0 * q_factor)

b0 =  +1.0 - alpha
b1 =  -2.0 * np.cos(w0)
b2 =  +1.0 + alpha
a0 =  +1.0 + alpha
a1 =  -2.0 * np.cos(w0)
a2 =  +1.0 - alpha

return (b0,b1,b2,a0,a1,a2)

# Constant 0 dB peak gain
# FIXME: Results in a peaking filter at freq1 rather than BP filter that's flat within the pass-band
def _make_biquad_bandpass( filter_freq, filter_freq2 ):

w0 = 2.0 * np.pi * filter_freq
BW = filter_freq2 - filter_freq
alpha = np.sin(w0) * np.sinh( np.log(2)/2 * BW * w0/np.sin(w0) )

b0 =  np.sin(w0) / 2.0
b1 =  +0.0
b2 = -np.sin(w0) / 2.0
a0 =  +1.0 + alpha
a1 =  -2.0 * np.cos(w0)
a2 =  +1.0 - alpha

return (b0,b1,b2,a0,a1,a2)

# gain can be + or -
def _make_biquad_peaking( filter_freq, q_factor, gain_db ):

A  = np.sqrt( 10 ** (gain_db/20) )
w0 = 2.0 * np.pi * filter_freq
alpha = np.sin(w0)/(2.0 * q_factor)

b0 = +1.0 + alpha * A
b1 = -2.0 * np.cos(w0)
b2 = +1.0 - alpha * A
a0 = +1.0 + alpha / A
a1 = -2.0 * np.cos(w0)
a2 = +1.0 - alpha / A

return (b0,b1,b2,a0,a1,a2)

def _make_biquad_lowshelf( filter_freq, q_factor, gain_db ):

S = q_factor
A  = 10.0 ** (gain_db / 40.0)
w0 = 2.0 * np.pi * filter_freq
alpha = np.sin(w0)/2 * np.sqrt( (A + 1/A)*(1/S - 1) + 2 )

b0 =    A*( (A+1) - (A-1)*np.cos(w0) + 2*np.sqrt(A)*alpha )
b1 =  2*A*( (A-1) - (A+1)*np.cos(w0)                   )
b2 =    A*( (A+1) - (A-1)*np.cos(w0) - 2*np.sqrt(A)*alpha )
a0 =        (A+1) + (A-1)*np.cos(w0) + 2*np.sqrt(A)*alpha
a1 =   -2*( (A-1) + (A+1)*np.cos(w0)                   )
a2 =        (A+1) + (A-1)*np.cos(w0) - 2*np.sqrt(A)*alpha

return (b0,b1,b2,a0,a1,a2)

def _make_biquad_highshelf( filter_freq, q_factor, gain_db ):

S = q_factor
A  = 10.0 ** (gain_db / 40.0)
w0 = 2.0 * np.pi * filter_freq
alpha = np.sin(w0)/2 * np.sqrt( (A + 1/A)*(1/S - 1) + 2 )

b0 =    A*( (A+1) + (A-1)*np.cos(w0) + 2*np.sqrt(A)*alpha )
b1 = -2*A*( (A-1) + (A+1)*np.cos(w0)                   )
b2 =    A*( (A+1) + (A-1)*np.cos(w0) - 2*np.sqrt(A)*alpha )
a0 =        (A+1) - (A-1)*np.cos(w0) + 2*np.sqrt(A)*alpha
a1 =    2*( (A-1) - (A+1)*np.cos(w0)                   )
a2 =        (A+1) - (A-1)*np.cos(w0) - 2*np.sqrt(A)*alpha

return (b0,b1,b2,a0,a1,a2)

def make_fir( type, frequency, order, window ):

taps = order + 1
#bb = signal.firwin( taps, cutoff=frequency, window="hamming" )
bb = dsp.firwin( taps, cutoff=frequency, window=window )
return (bb,[1])

def make_iir( frequency, type, order, shape ):

#for i in range(0,len(fc)): fc[i] = 2 * fc[i] / fs
return signal.iirfilter( order, frequency, btype=type, ftype=shape )

def make_biq( type, frequency, q_factor, gain_db=1.0 ):

freq = float( frequency )
Q    = float( q_factor )
gain = float( gain_db )

if type == "notch":      (b0,b1,b2,a0,a1,a2) = _make_biquad_notch    ( freq, Q )
if type == "lowpass":    (b0,b1,b2,a0,a1,a2) = _make_biquad_lowpass  ( freq, Q )
if type == "highpass":   (b0,b1,b2,a0,a1,a2) = _make_biquad_highpass ( freq, Q )
if type == "allpass":    (b0,b1,b2,a0,a1,a2) = _make_biquad_allpass  ( freq, Q )
if type == "bandpass":   (b0,b1,b2,a0,a1,a2) = _make_biquad_bandpass ( freq, Q )
if type == "peaking":    (b0,b1,b2,a0,a1,a2) = _make_biquad_peaking  ( freq, Q, gain )
if type == "highshelf":  (b0,b1,b2,a0,a1,a2) = _make_biquad_highshelf( freq, Q, gain )
if type == "lowshelf":   (b0,b1,b2,a0,a1,a2) = _make_biquad_lowshelf ( freq, Q, gain )

return ([b0,b1,b2],[a0,a1,a2])

if __name__ == "__main__":

    _gain = 1.0
    _fmt = ""
    _plot = False
    _xmin = None
    _xmax = None
    _ymin = None
    _ymax = None
    _window = "nuttall" # "hamming"
    _name = None

    bb = None
    aa = None

    for i in range( 1, len(sys.argv) ):

        option = sys.argv[i].lower()

        if option.find("name=") == 0:   _name  = option[option.find("name=")+len("name="):]
        if option.find("freq=") == 0:   _freq  = float( option[option.find("freq=")+len("freq="):] )
        if option.find("freq2=") == 0:  _Q     = float( option[option.find("freq2=")+len("freq2="):] )
        if option.find("type=") == 0:   _type  = option[option.find("type=")+len("type="):]
        if option.find("shape=") == 0:  _shape = option[option.find("shape=")+len("shape="):]
        if option.find("gain=") == 0:   _gain  = float( option[option.find("gain=")+len("gain="):] )
        if option.find("q=") == 0:      _Q     = float( option[option.find("q=")+len("q="):] )
        if option.find("order=") == 0:  _order = int( option[option.find("order=")+len("order="):] )
        if option.find("format=") == 0: _fmt   = option[option.find("format=")+len("format="):]
        if option.find("plot=") == 0:   _plot  = option[option.find("plot=")+len("plot="):] == "y"
        if option.find("plot=") == 0:   _plot  = option[option.find("plot=")+len("plot="):] == "yes"
        if option.find("xmin=") == 0:   _xmin  = float( option[option.find("xmin=")+len("xmin="):] )
        if option.find("xmax=") == 0:   _xmax  = float( option[option.find("xmax=")+len("xmax="):] )
        if option.find("ymin=") == 0:   _ymin  = float( option[option.find("ymin=")+len("ymin="):] )
        if option.find("ymax=") == 0:   _ymax  = float( option[option.find("ymax=")+len("ymax="):] )

    if sys.argv[1].lower() == "fir":

        options = "FIR Type=\"%s\" Fc=%u Order=\"%u\" Window=\"%s\"" % (_type,_freq,_order,_window)
        (bb,aa) = make_fir( _type, _freq, _order, _window )
        if _plot: plot_response( bb, aa, _xmin, _xmax, _ymin, _ymax )
        options += " ---> b0,...b%u" % (len(bb)-1)

    if sys.argv[1].lower() == "iir":

        options = "IIR Type=\"%s\" Fs=%u Fc=%u Order=\"%u\" Shape=\"%s\"" % (_type,_freq,_order,_shape)
        (bb,aa) = make_iir( _fc / (_fs/2.0), _type, _order, _shape )
        if _plot: plot_response( bb, aa, _fs / 2.0, _xmin, _xmax, _ymin, _ymax )
        options += " ---> b0,...b%u,a1,...a%u" % (len(bb)-1,len(aa)-1)

    if sys.argv[1].lower() == "biquad":

        options = "BiQuad Type=\"%s\" Freq=%f Q=%04.3f Gain=%03.1f" % (_type,_freq,_Q,_gain)
        (bb,aa) = make_biq( _type, _freq, _Q, _gain )
        if _plot: plot_response( bb, aa, _xmin, _xmax, _ymin, _ymax )
        options += " ---> b0,...b%u,a1,...a%u" % (len(bb)-1,len(aa)-1)

    if _name == None and _fmt != None:
        print options
        print coefficients_string( _fmt.upper(), bb, aa )


And finally, a plot of the filter magnitude and phase response ...



markseel

If you want to do some non-real time DSP here's some more scripts  :D

One script parses a WAVE file, extracts the data and prints the values to STDOUT.  Multiple WAVE files are supported - the output has as many columns as WAVE files specified for parsing.  The other script plots data in the time or frequency domain.

Example:  Parse two impulse response WAVE files and direct output to a data file.  Then plot the data.  This plot shows the frequency response of the two IR sample files that were just parsed ...


python wave2data.py ir1.wav ir2.wav > data
python plot.py data freq log


Here's the two IR responses.  One is for a 1x12 speaker/cab and the other is for a 4x12 speaker cab with some room dynamics included ... very different responses!



WAVE parsing script ...


import sys, struct

def _assert( condition, message ):

    if not condition:
        print message
        exit( 0 )

def _parse_wave( wave_file ):

    rate = 0; samples = None
    (group_id,total_size,type_id) = struct.unpack( "<III", wave_file.read(12) )
    _assert( group_id == 0x46464952, "Unknown File Format" ) # Signature for 'RIFF'
    _assert( type_id  == 0x45564157, "Unknown File Format" ) # Signature for 'WAVE'
    while True:
        if total_size <= 8: break
        data = wave_file.read(8)
        if len(data) < 8: break;
        (blockid,blocksz) = struct.unpack( "<II", data )
        #print "WaveIn: BlockID=0x%04X BlockSize=%u" % (blockid,blocksz)
        #print "WaveIn: ByteCount=%u" % (blocksz)
        total_size -= 8
        if blockid == 0x20746D66: # Signature for 'fmt'
            if total_size <= 16: break
            (format,channels,rate,thruput,align,width) = struct.unpack( "<HHIIHH", wave_file.read(16) )
            #print "WaveIn: ByteCount=%u Channels=%u Rate=%u WordSize=%u" % (blocksz,channels,rate,width)
            #print "WaveIn: Format=%u Alignment=%u" % (format,align)
            total_size -= 16
        elif blockid == 0x61746164: # Signature for 'data'
            samples = [0] * (blocksz / (width/8))
            count = 0
            data = wave_file.read( blocksz )
            if channels == 1:
                while len(data) >= width/8:
                    if width == 8:
                        samples[count] = struct.unpack( "b", data[0:1] )[0] * 255 * 256 * 256
                        data = data[1:]
                    if width == 16:
                        samples[count]  = struct.unpack( "b", data[1:2] )[0] * 255 * 256 * 256
                        samples[count] += struct.unpack( "B", data[0:1] )[0] * 255 * 256
                        data = data[2:]
                    if width == 24:
                        samples[count]  = struct.unpack( "b", data[2:3] )[0] * 255 * 256 * 256
                        samples[count] += struct.unpack( "B", data[1:2] )[0] * 256 * 256
                        samples[count] += struct.unpack( "B", data[0:1] )[0] * 256
                        data = data[3:]
                    if width == 32:
                        samples[count]  = struct.unpack( "b", data[3:4] )[0] * 255 * 256 * 256
                        samples[count] += struct.unpack( "B", data[2:3] )[0] * 256 * 256
                        samples[count] += struct.unpack( "B", data[1:2] )[0] * 256
                        samples[count] += struct.unpack( "B", data[0:1] )[0]
                        data = data[4:]
                    count += 1
            total_size -= blocksz
    return (rate,samples)

data = []
count = 99999999

for name in sys.argv[1:]:

    try: file = open( name, "rb" )
    except IOError as err: file = None
    _assert( file != None, "Unable to open file" )
    if file != None:
        (rate,samples) = _parse_wave( file )
        if count > len(samples): count = len(samples)
        data.append( samples )
        file.close()

for row in range(0,count):
    for col in range(0,len(data)):
        val = data[col][row]
        sys.stdout.write( "%08x " % (data[col][row] & 0xFFFFFFFF) )
    sys.stdout.write( "\n" )


Data plotting script ...


import sys, struct, random

print ""
print "Usage: python plot.py <datafile> time"
print "       python plot.py <datafile> time [beg end]"
print "       python plot.py <datafile> freq lin"
print "       python plot.py <datafile> freq log"
print ""
print "Where: <datafile> Contains one sample value per line.  Each sample is an"
print "                  ASCII/HEX value (e.g. FFFF0001) representing a fixed-"
print "                  point sample value."
print "       time       Indicates that a time-domain plot should be shown"
print "       freq       Indicates that a frequency-domain plot should be shown"
print "       [beg end]  Optional; specifies the first and last sample in order to"
print "                  create a sub-set of data to be plotted"
print ""
print "Examle: Create time-domain plot data in \'out.txt\' showing samples 100"
print "        through 300 ... bash$ python plot.py out.txt 100 300"
print ""
print "Examle: Create frequency-domain plot data in \'out.txt\' showing the Y-axis"
print "        with a logarithmic scale ... bash$ python plot.py out.txt freq log"
print ""

if len(sys.argv) < 3: exit(0)

import numpy
import matplotlib.pyplot

numpy.seterr( all='ignore' )

K = 1
yy = [[0] * 65536, [0] * 65536, [0] * 65536]
mm = [[0] * 65536, [0] * 65536, [0] * 65536]
N = 0

file = open( sys.argv[1], "rt" )
while True:
    ll = file.readline();
    ll = ll.replace( "\n", "" );
    if len(ll) < 1: break
    ll = ll.split()
    K = len(ll)
    for ii in range(0,K):
        xx = int(ll[ii],16) * 2
        tt = xx
        if xx & 0x80000000 == 0x80000000:
            xx = -(1 - float((xx & 0x7FFFFFFF)) / (2**31))
        else:
            xx = float((xx & 0x7FFFFFFF)) / ((2**31)-1)
        yy[ii][N] = xx
    N += 1
file.close()

xx = [0] * N
ff = [0] * N
mm[0] = [0] * N
mm[1] = [0] * N
mm[2] = [0] * N
yy[0] = yy[0][0:N]
yy[1] = yy[1][0:N]
yy[2] = yy[2][0:N]

for ii in range( 0, N ):
    ff[ii] = ii
    xx[ii] = float(ii) / N

matplotlib.pyplot.grid( b=True, which='both', color='0.65',linestyle='-' )

if sys.argv[2] == "time":

    beg = 0; end = N
    if len(sys.argv) > 3: beg = int( sys.argv[3] )
    if len(sys.argv) > 4: end = int( sys.argv[4] )
    ff = ff[beg:end]

    matplotlib.pyplot.title( "Time Domain" )
    matplotlib.pyplot.xlabel( "Sample Number" )
    matplotlib.pyplot.ylabel( "Sample Amplitude" )
    matplotlib.pyplot.xticks( numpy.arange(min(ff), max(ff), (max(ff)-min(ff))/4 ))

    ymin = 0; ymax = 0

    for ii in range(0,K):
        yy[ii] = yy[ii][beg:end]
        ymin = min(ymin,min(yy[ii]))
        ymax = max(ymax,max(yy[ii]))
        matplotlib.pyplot.plot( ff, yy[ii] )

    matplotlib.pyplot.ylim( ymin, ymax )
    matplotlib.pyplot.yticks( numpy.arange(ymin, ymax, (ymax-ymin)/20 ))
    matplotlib.pyplot.show()

if sys.argv[2] == "freq":

    fs = 1.0
    if len(sys.argv) > 4: fs = float( sys.argv[4] )

    ff = numpy.array(ff[0:N/2+1]) / float(N/2)
    ff *= fs / 2.0

    for ii in range(0,K):
        mm[ii] = numpy.fft.rfft(yy[ii])
        mm[ii] = mm[ii][0:len(mm[ii])-0]
        #mm[ii] = numpy.abs(mm[ii]) * 2.0 / N
        mm[ii] = numpy.abs(mm[ii]) * 1.0 / N
        mm[ii] /= numpy.max(mm[ii])
        mm[ii] = mm[ii][0:N/2+1]

    mode = "lin"
    if len(sys.argv) > 3:
        if( sys.argv[3] == "log" ): mode = "log"

    if mode == "lin":

        matplotlib.pyplot.title( "Frequency Domain" )
        matplotlib.pyplot.xlabel( "Frequency (Normalized to Fs)" )
        matplotlib.pyplot.ylabel( "FFT Amplitude (Normalized to Max)" )
        matplotlib.pyplot.ylim( 0, max(mm) )
        matplotlib.pyplot.plot( ff, mm )
        matplotlib.pyplot.xticks( numpy.arange(min(ff), max(ff), (max(ff)-min(ff))/10 ))
        matplotlib.pyplot.yticks( numpy.arange(0, max(mm), max(mm)/10.0 ))
        matplotlib.pyplot.show()

    if mode == "log":

        matplotlib.pyplot.title( "Frequency Domain" )
        matplotlib.pyplot.xlabel( "Frequency (Normalized to Fs)" )
        matplotlib.pyplot.ylabel( "FFT Amplitude (dBfs)" )
        matplotlib.pyplot.ylim( -130, 20 )

        for ii in range(0,K):
            mm[ii] = 20.0 * numpy.log( mm[ii] )
            matplotlib.pyplot.plot( ff, mm[ii] )

        matplotlib.pyplot.xticks( numpy.arange(min(ff), max(ff), (max(ff)-min(ff))/10 ))
        matplotlib.pyplot.yticks( numpy.arange(-130, 20, (20+130)/15 ))
        matplotlib.pyplot.show()


potul

if you are willing to program in C, a good option is the project

http://www.diystompboxes.com/smfforum/index.php?topic=116415.0

As it uses a dsPIC, it can be programmed with C.

Mat

will20fitz

Thanks all for your responses...

Plenty of food for thought... but its surprising how quickly you feel out of your depth with DSP! :-)

Thats never stopped me before though!

Let me do some more research and i will come back with some questions!

Cheers again everyone...