Mailman 3 February 2004 - SciPy-User

help and lazy-import (?)
by Arnd Baecker March 2, 2004

March 2, 2004

Hi, there there seems to be something weird going on when trying to access help in the following way >>> import scipy >>> help("scipy.linalg") Help on instance of _ModuleLoader in scipy: scipy.linalg = <module 'scipy.linalg' from '/opt/python/lib/pyt... e-packages/scipy/linalg/__init__.pyc' [imported]> Doing another >>> help("scipy.linalg") one gets the expected help. I'd guess that this somehow has to to with the lazy importer of scipy. But maybe there is a workaround? Many thanks in advance, Arnd

3 6

RE: [SciPy-user] FIR Filters
by Bob.Cowdery＠CGI-Europe.com Feb. 27, 2004

Feb. 27, 2004

Travis With your help I am making progress. I started off using numarray, then ran into problems so I am now using Numeric. I only import scipy so I am using whatever it imports which I believe is Numeric. I have tried to follow your advice and move to complex types as I can see it makes sense. I don't think I have done a complete job yet. The audio output has a pulse type modulation on it which should not be there, however selecting different passbands is definitely having the desired effect. I have copied all the code in below for context. Yes the filtering is in the frequency domain which I am told gives better stopband attenuation and good shape. The way the filter and demodulation works is still a bit of a mystery to me as I didn't write the original code. I would love to work out exactly what it does but I'm not quite there yet. I can't really ask very targetted questions yet until I understand more. The code cleary still has problems but I don't know why or where. The good news is the processor is down to 50% now from 80%. If you do have a few moments perhaps you could point out any more obvious errors. Thanks Bob import time import threading import fastaudio as f from scipy import * class dsp( threading.Thread ): # signal processing def sig( self, samples ) : # allocate the filter arrays global m_filterMagn global m_filterPhase # convert to a numeric array of type double x = asarray(samples, typecode=Float64) # get a complex array ordered from the I,Q source complexin = x[1::2] + 1j*x[::2] # convert to complex frequency domain cart_complex_in_fd = fft(complexin) # apply filter to magnitude and phase magn = abs(complexin)*m_filterMagn phase = angle(complexin)*m_filterPhase # polar to cartesian in frequency domain demodulated = magn * cos(phase) + 1j*(magn * sin(phase)) # convert to time domain complexout = ifft(demodulated) # convert back to interleaved array x[1::2] = real(complexout) x[::2] = real(complexout) # convert back to an integer array x = x.astype('l') # return sequence of processed samples return x def firwin(self, N, cutoff, width=None, window='blackman'): """ FIR Filter Design using windowed ideal filter method. Inputs: N -- order of filter (number of taps) cutoff -- cutoff frequency of filter (normalized so that 1 corresponds to Nyquist or pi radians / sample) width -- if width is not None, then assume it is the approximate width of the transition region (normalized so that 1 corresonds to pi) for use in kaiser FIR filter design. window -- desired window to use. Outputs: h -- coefficients of length N fir filter. """ if isinstance(width,float): A = 2.285*N*width + 8 if (A < 21): beta = 0.0 elif (A <= 50): beta = 0.5842*(A-21)**0.4 + 0.07886*(A-21) else: beta = 0.1102*(A-8.7) window=('kaiser',beta) win = signal.get_window(window,N,fftbins=1) alpha = N//2 m = signal.Num.arange(0,N) return win*special.sinc(cutoff*(m-alpha)) # calculate the complex filter def filter_calc( self, fLow, fHigh ): global m_filterMagn global m_filterPhase m_filterMagn = zeros(2048) m_filterPhase = zeros(2048) # create empty arrays to hold results Rh = zeros(2048) Ih = zeros(2048) # calculate the high/low cutoff freq as a fraction of the sample rate fhc = fHigh/44100.0 flc = fLow/44100.0 # calculate fft bin size bin_size = 44100.0/4096.0 # calculate no of bins in filter width bins = fHigh/bin_size # calculate FIR lowpass filter coefficients Rh = self.firwin(2048.0,(fhc-flc)/2.0) # compute USB filter by modulating the lowpass filter up filterPhase = 0.0 filterFreq = (flc+fhc)*pi phaseSin = 0 phaseCos = 0 # modulate to USB centre frequency i = 0 while i< 2048 : phaseSin = sin(filterPhase) phaseCos = cos(filterPhase) Ih[i] = Rh[i]*phaseSin Rh[i] = Rh[i]*phaseCos filterPhase = filterPhase+filterFreq if filterPhase >= 6.28318530717959 : filterPhase = filterPhase - 6.28318530717959 i = i+1 # for LSB create conjugate of USB filter # comment out for LSB Ih = -Ih # compute complex frequency domain response of LSB or USB filter complexin = Rh + 1j*Ih complexfd = fft(complexin) # Cartesian to polar gives the magnitude and phase components of the fft result m_filterMagn = abs(complexfd) m_filterPhase = angle(complexfd) # start streaming def run( self ): global running running = 1 #s = f.stream(7500, # sample rate # 2, # number of channels # 'int16', # sample format, choose from 'int8', 'int16', 'int32'. # 4096, # frames per buffer # 16) # maximum number of input buffers # open stream for input/output print "Create, open and start stream" global s s = f.stream(44100,2,'int16',2048,4) # open and start streaming s.open() s.start() print "Thread running..." # calculate the SSB filter self.filter_calc(300.0,3000.0) while( running ): # read next block of data from the input buffer samples = s.read() if( len( samples ) > 0 ): demod = self.sig( samples ) # now, write the data back s.write( demod ) # sleep for 10 ms # a stero block of 4096 samples takes 46 ms so 10ms should avoid dropped blocks time.sleep(0.01) print "Thread exiting..." def start(): global DSP DSP = dsp() DSP.start() def stop(): # stop the stream global running running = 0 global s s.stop() # and close it s.close() def filter(id): # set the filter passband global DSP if id == 0: DSP.filter_calc(300.0,3000.0) elif id == 1: DSP.filter_calc(300.0,2000.0) elif id == 2: DSP.filter_calc(500.0,700.0) -----Original Message----- From: Travis Oliphant [mailto:oliphant@ee.byu.edu] Sent: 26 February 2004 22:38 To: SciPy Users List Subject: Re: [SciPy-user] FIR Filters Bob.Cowdery(a)CGI-Europe.com wrote: > Travis > > Thanks very much for the responses. Don't confuse me with someone who > knows what they are doing (but I am working on it!). To give you some > background I am converting a .NET application that uses Intel IPP to > Python. The application is for a software defined radio. At the moment > I have audio capture which contains my quadrature samples in blocks of > 4096 (2048 each for I&Q) with 44100 sampling rate (46ms per block). I > have a simple dsp process to test each signal processing call. At the > moment this should output exactly what it receives, i.e just the audio > input as the demodulator needs the FIR filter and some magic applied > at the point indicated in the code. When I have written the > demodulator I can test firwin() properly. At the moment I can see it > outputs valid data, and it's fast which is good. The code below is the > very simple path. I have had to write two functions for > polar/cartesian conversions as I didn't find any in the library. I > don't know if these work right as I can hear clicking sounds on the > audio not there when they are commented out. However, it could be > audio overrun as I am up at 75% processor. The functions are also > quite slow (adding 30% to > > my processor loading). Is there a more efficient way to do this? > All the best in your education... First, Are you using numarray? if so then that might cause a slower application than using Numeric (older but often faster arraymodule) Comments on your code below: > def sigproc( self, samples ) : > > # convert to a numarray array of type double > x = asarray(samples, typecode=Float64) > is this numarray.asarray or Numeric.asarray ? > > # first to do is sort into real and image arrays from > interleaved array > realin = x[1::2] # every other element starting at index 1 > imagin = x[::2] # every other element starting at index 0 \ > > > # convert to frequency domain > realfd = fft(realin) > imagfd = fft(imagin) > Normally, real and imaginary parts of the data are handled together when a complex data type is available. I would do in = x[1::2] + 1j*x[::2] fd = fft(in) > > # cartesian to polar in frequency domain > magn,phase = self.cartesian_to_polar(realfd, imagfd) > > This seems strange to me, why are you filtering the magnitude and phase separately? This is very non-standard. > > # apply filter and demodulate > I doubt you want a 2048 tap filter at this point. Probably something a lot less. What kind of modulation are you demodulating? Are you doing frequency-domain filtering then? > > # polar to cartesian in frequency domain > > realfd, imagfd = self.polar_to_cartesian(magn, phase) > > > # convert to time domain > realtd = ifft(realfd) > imagtd = ifft(imagfd) > > # convert back to interleaved array > x[::2] = realtd.astype(Float64) > x[1::2] = imagtd.astype(Float64) > > # convert back to an integer array > x = x.astype('l') > > # return python list of processed samples > return x.tolist() > Why are you going back to a python list? A numeric array is a valid sequence. Why does it need to be a list. You have a lot of memory copying going on which is going to slow down everything. > > # utility functions not in the library > def polar_to_cartesian(self, magn, phase): > return magn * cos(phase), magn * sin(phase) > > def cartesian_to_polar(self, real, imag): > # convert the arrays real, imag into the polar system > # magnitude array > magn = sqrt(real**2 + imag**2) > # phase array > phase = > choose(equal(real,0),(arctan(imag/real),(choose(less(imag,0),((pi/2),-(pi/2) ))))) > > return magn, phase > > If you have a complex-valued array, then abs(x) is the magnitude and angle(x) is the phase. My colleagues down the hall do software radio and I have a lot of background in signal processing. Let me know if you have more specific questions. Thanks, -Travis O. _______________________________________________ SciPy-user mailing list SciPy-user(a)scipy.net http://www.scipy.net/mailman/listinfo/scipy-user *** Confidentiality Notice *** Proprietary/Confidential Information belonging to CGI Group Inc. and its affiliates may be contained in this message. If you are not a recipient indicated or intended in this message (or responsible for delivery of this message to such person), or you think for any reason that this message may have been addressed to you in error, you may not use or copy or deliver this message to anyone else. In such case, you should destroy this message and are asked to notify the sender by reply email.

1 0

Boa Constructor Plug-in for ploting
by Ricardo Henriques Feb. 27, 2004

Feb. 27, 2004

Is there a Boa Constructor plug-in using a scipy based plot interface like Chaco or scipy.plt out there? I have been using wxPyPlot witch is based on scipy.plot but laks many of its features... Tks... Paxcal

2 1

RE: [SciPy-user] FIR Filters
by Bob.Cowdery＠CGI-Europe.com Feb. 26, 2004

Feb. 26, 2004

Travis Thanks very much for the responses. Don't confuse me with someone who knows what they are doing (but I am working on it!). To give you some background I am converting a .NET application that uses Intel IPP to Python. The application is for a software defined radio. At the moment I have audio capture which contains my quadrature samples in blocks of 4096 (2048 each for I&Q) with 44100 sampling rate (46ms per block). I have a simple dsp process to test each signal processing call. At the moment this should output exactly what it receives, i.e just the audio input as the demodulator needs the FIR filter and some magic applied at the point indicated in the code. When I have written the demodulator I can test firwin() properly. At the moment I can see it outputs valid data, and it's fast which is good. The code below is the very simple path. I have had to write two functions for polar/cartesian conversions as I didn't find any in the library. I don't know if these work right as I can hear clicking sounds on the audio not there when they are commented out. However, it could be audio overrun as I am up at 75% processor. The functions are also quite slow (adding 30% to my processor loading). Is there a more efficient way to do this? def sigproc( self, samples ) : # convert to a numarray array of type double x = asarray(samples, typecode=Float64) # first to do is sort into real and image arrays from interleaved array realin = x[1::2] # every other element starting at index 1 imagin = x[::2] # every other element starting at index 0 # convert to frequency domain realfd = fft(realin) imagfd = fft(imagin) # cartesian to polar in frequency domain magn,phase = self.cartesian_to_polar(realfd, imagfd) # apply filter and demodulate # polar to cartesian in frequency domain realfd, imagfd = self.polar_to_cartesian(magn, phase) # convert to time domain realtd = ifft(realfd) imagtd = ifft(imagfd) # convert back to interleaved array x[::2] = realtd.astype(Float64) x[1::2] = imagtd.astype(Float64) # convert back to an integer array x = x.astype('l') # return python list of processed samples return x.tolist() # utility functions not in the library def polar_to_cartesian(self, magn, phase): return magn * cos(phase), magn * sin(phase) def cartesian_to_polar(self, real, imag): # convert the arrays real, imag into the polar system # magnitude array magn = sqrt(real**2 + imag**2) # phase array phase = choose(equal(real,0),(arctan(imag/real),(choose(less(imag,0),((pi/2),-(pi/2) ))))) return magn, phase -----Original Message----- From: Travis Oliphant [mailto:oliphant@ee.byu.edu] Sent: 26 February 2004 07:10 To: SciPy Users List Subject: Re: [SciPy-user] FIR Filters > I am trying to get a windowed FIR filter to use in fast convolution > filtering. > > I have tried using : > signal.remez(2048,[300,3000],[1],Hz=44100) > to generate a 2048 tap filter, which is what I have used previously in > Intel IPP. Firstly, this gives me an array full of QNAN errors. Infact > the only time I don't get errors is by using a very small number of > taps, say 10. Various other values of taps fails to converge or gives > QNAN errors. With the 2048 taps it is also very, very slow (about 75 > seconds on a 1.2GHz). The IPP implementation is sub-second. I don't > know if it's throwing an exception on every QNAN or quite what is > hapenning. > > I am also confused by not being able to specify a window type and I > really need lowpass not bandpass. I can see there is a function > get_window() but I don't know how this relates to an FIR filter when > the only other parameter is fftbins. Do I do something like tell it > the number of fftbins I want the filter to cover (I have 2048 bins of > around 10Hz each) then shift it to the correct part of the spectrum? Bob, Here is a simple windowed ideal lowpass filter algorithm. The two important parameters are the number of taps and the cutoff frequency normalized so that 1 corresponds to pi radians/sample (i.e. the Nyquist frequency). If you specify width it will choose a nice Kaiser window for you to try and reach that width of transistion region (from passband to stopband). Otherwise, you can specify the window type (see signal.get_window for types) This is now in SciPy CVS. def firwin(N, cutoff, width=None, window='hamming'): """FIR Filter Design using windowed ideal filter method. Inputs: N -- order of filter (number of taps) cutoff -- cutoff frequency of filter (normalized so that 1 corresponds to Nyquist or pi radians / sample) width -- if width is not None, then assume it is the approximate width of the transition region (normalized so that 1 corresonds to pi) for use in kaiser FIR filter design. window -- desired window to use. Outputs: h -- coefficients of length N fir filter. """ from signaltools import get_window if isinstance(width,float): A = 2.285*N*width + 8 if (A < 21): beta = 0.0 elif (A <= 50): beta = 0.5842*(A-21)**0.4 + 0.07886*(A-21) else: beta = 0.1102*(A-8.7) window=('kaiser',beta) win = get_window(window,N,fftbins=1) alpha = N//2 m = Num.arange(0,N) return win*special.sinc(cutoff*(m-alpha)) _______________________________________________ SciPy-user mailing list SciPy-user(a)scipy.net http://www.scipy.net/mailman/listinfo/scipy-user *** Confidentiality Notice *** Proprietary/Confidential Information belonging to CGI Group Inc. and its affiliates may be contained in this message. If you are not a recipient indicated or intended in this message (or responsible for delivery of this message to such person), or you think for any reason that this message may have been addressed to you in error, you may not use or copy or deliver this message to anyone else. In such case, you should destroy this message and are asked to notify the sender by reply email.

2 1

[SciPy-user] Re: data smoothing: interpolate.splrep ignores s parameter
by george young Feb. 26, 2004

Feb. 26, 2004

[[reposting with "rough" data file appended, sorry]] [SciPy-0.2.0_alpha_200.4161, Numeric-23.1, Python 2.3.3, x86 Linux] My goal really is to smooth some noisy measurement data without messing up it's *shape*. My first attempt was 1d splines. I did: from Numeric import * import scipy f = file('rough') data = array([(float(x),float(y)) for (x,y) in [l.split() for l in f]]) rep = scipy.interpolate.splrep(data[:,0], data[:,1], s=s) smooth_y = scipy.interpolate.splev(data[:,0], rep) "s" is supposed to vary the amount of smoothing. For s=0, i get my original data back, as expected. But for all other values, including the recommended range of (m-sqrt(2.0*m)) to (m+sqrt(2.0*m), I get a single, too much smoothed, result. It seems like splrep is not using the "s" value to adjust the splines, except to sense that it's not zero. Splines may not be the right method anyway, since they tend to warp the shape of the curve, and I need to get the data's derivatives. Is there some way to fashion a low pass filter? It seems like fft should be useful here, but I have very little experience with fft's. -- George Young "rough" file: 0.5 -1.7355e-10 0.455 -4.6185e-10 0.41 -3.2125e-10 0.365 -4.035e-10 0.32 2.951e-10 0.275 2.628e-10 0.23 2.922e-10 0.185 -3.1775e-10 0.14 -3.3209e-10 0.095 -3.6665e-10 0.05 -2.7394e-10 0.0 2.98e-10 -0.04 3.228e-10 -0.085 3.219e-10 -0.13 4.093e-10 -0.175 2.7406e-10 -0.22 1.102e-10 -0.265 -6.012e-10 -0.31 -3.9674e-09 -0.355 -1.2145e-08 -0.4 -3.9299e-08 -0.445 -1.2074e-07 -0.49 -3.2109e-07 -0.535 -7.3798e-07 -0.58 -1.4695e-06 -0.625 -2.548e-06 -0.67 -3.9104e-06 -0.715 -5.4895e-06 -0.76 -7.1758e-06 -0.805 -8.9267e-06 -0.85 -1.0675e-05 -0.895 -1.2425e-05 -0.94 -1.4115e-05 -0.985 -1.579e-05 -1.03 -1.7415e-05 -1.075 -1.9005e-05 -1.12 -2.0505e-05 -1.165 -2.2e-05 -1.21 -2.346e-05 -1.255 -2.4885e-05 -1.3 -2.624e-05 -1.345 -2.757e-05 -1.39 -2.886e-05 -1.435 -3.015e-05 -1.48 -3.1404e-05 -1.525 -3.2604e-05 -1.57 -3.385e-05 -1.615 -3.498e-05 -1.66 -3.6105e-05 -1.705 -3.718e-05 -1.75 -3.8279e-05 -1.795 -3.932e-05 -1.84 -4.033e-05 -1.885 -4.13e-05 -1.93 -4.2304e-05 -1.975 -4.3275e-05 -2.02 -4.429e-05 -2.065 -4.5205e-05 -2.11 -4.6069e-05 -2.155 -4.697e-05 -2.2 -4.7825e-05 -2.245 -4.8639e-05 -2.29 -4.9515e-05 -2.335 -5.0319e-05 -2.38 -5.1094e-05 -2.425 -5.188e-05 -2.47 -5.264e-05 -2.515 -5.3454e-05 -2.56 -5.4169e-05 -2.605 -5.4905e-05 -2.65 -5.5633e-05 -2.695 -5.6334e-05 -2.74 -5.7043e-05 -2.785 -5.7749e-05 -2.83 -5.8343e-05 -2.875 -5.9085e-05 -2.92 -5.9729e-05 -2.965 -6.0344e-05 -3.01 -6.094e-05 -3.055 -6.1564e-05 -3.1 -6.2104e-05 -3.145 -6.2749e-05 -3.19 -6.3348e-05 -3.235 -6.3956e-05 -3.28 -6.4414e-05 -3.325 -6.4958e-05 -3.37 -6.555e-05 -3.415 -6.6075e-05 -3.46 -6.656e-05 -3.505 -6.71e-05 -3.55 -6.7603e-05 -3.595 -6.8083e-05 -3.64 -6.8623e-05 -3.685 -6.9048e-05 -3.73 -6.9518e-05 -3.775 -6.998e-05 -3.82 -7.0508e-05 -3.865 -7.0929e-05 -3.91 -7.1343e-05 -3.955 -7.1824e-05 -4.0 -7.2218e-05

5 5

scipy cvs snapshot problem
by Jim Benson Feb. 26, 2004

Feb. 26, 2004

Hi, I downloaded a cvs snapshot on 2004-02-25 and installed it over a previous cvs snapshot (from 2003-12-03). The previous one did pass all scipy.test levels. The recent version: Found scipy_distutils version '0.2.2_alpha_26.374' in '/usr/local/lib/python2.3/ site-packages/scipy_distutils/__init__.pyc' Importing scipy_distutils.command.build_flib ... ok failed the first test: ====================================================================== FAIL: check_nils_20Feb04 (scipy.linalg.basic.test_basic.test_solve) ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/local/lib/python2.3/site-packages/scipy/linalg/tests/test_basic.py", line 112, in check_nils_20Feb04 assert_array_almost_equal(X,Ainv) File "/usr/local/lib/python2.3/site-packages/scipy_test/testing.py", line 636, in assert_array_almost_equal assert cond,\ AssertionError: Arrays are not almost equal (mismatch 50.0%): Array 1: [[-2.3261187+2.4103813j 2.6175944-0.5570154j] [ 2.9991264-1.7766662j -2.0232508-0.2587063j]] Array 2: [[-2.8831341-0.2072131j 2.6175944-0.5570154j] [ 2.7404201+0.2465847j -2.0232508-0.2587063j]] ---------------------------------------------------------------------- Ran 1198 tests in 1.918s FAILED (failures=1) <unittest.TextTestRunner object at 0x4240ef2c> >>> My atlas version: floyd:/usr/local/scipy/Lib/linalg>python -c 'import atlas_version' ATLAS version 3.4.2 built by root on Tue Dec 2 20:35:31 MST 2003: ...am i forgetting something critical? Thanks, Jim

3 2

weave, including headers in generated file
by Magne Westlie Feb. 26, 2004

Feb. 26, 2004

Hello. I'm having a problem with using weave which I was wondering if someone could help me out with. I've got a class Hello declared in a file hw.h, but how do I tell weave to include this file in the generated cpp-file? If I try: code = """ Hello h; h.hw(); """ weave.inline(code, [], compiler='gcc', verbose=2, force=1, sources=['/home/mw/weave/inline-ex/hw.cpp'], include_dirs=['/home/mw/weave/inline-ex/'], ) Compiling fails because Hello is not defined: /home/mw/.python23_compiled/sc_dbdc99b5e1f621cd206ee575625e18f15.cpp:644: error: ` Hello' undeclared (first use this function) If I open the cpp file manually and puts: #include <hw.h> in it, everything works. Magne Westlie

2 1

scipy.xplt.ghelp (2nd try ;-)
by Arnd Baecker Feb. 26, 2004

Feb. 26, 2004

Hi, this is a reformulation of my recent question with the same subject, I hope I get it right this time ;-) If I understand things correctly, scipy.xplt.ghelp is gone because the information from xplt/gistdata/help.help has been moved to the doc strings of the commands. If this is correct, how do I get information on the keywords for a plot command ? E.g. help("scipy.xplt.plg") works, but I don't know how to get information on the keywords, like `color`, `closed` etc. I think all the information is given here http://bonsai.ims.u-tokyo.ac.jp/~mdehoon/software/python/pygist_html/node57… How can I access this type of information from within python? (Before scipy.xplt.ghelp("color") did the job). Many thanks, Arnd

3 4

FIR Filters
by Bob.Cowdery＠CGI-Europe.com Feb. 26, 2004

Feb. 26, 2004

I am trying to get a windowed FIR filter to use in fast convolution filtering. I have tried using : signal.remez(2048,[300,3000],[1],Hz=44100) to generate a 2048 tap filter, which is what I have used previously in Intel IPP. Firstly, this gives me an array full of QNAN errors. Infact the only time I don't get errors is by using a very small number of taps, say 10. Various other values of taps fails to converge or gives QNAN errors. With the 2048 taps it is also very, very slow (about 75 seconds on a 1.2GHz). The IPP implementation is sub-second. I don't know if it's throwing an exception on every QNAN or quite what is hapenning. I am also confused by not being able to specify a window type and I really need lowpass not bandpass. I can see there is a function get_window() but I don't know how this relates to an FIR filter when the only other parameter is fftbins. Do I do something like tell it the number of fftbins I want the filter to cover (I have 2048 bins of around 10Hz each) then shift it to the correct part of the spectrum? Any help appreciated. Bob Bob Cowdery CGI Senior Technical Architect +44(0)1438 791517 Mobile: +44(0)7771 532138 bob.cowdery(a)cgi-europe.com <mailto:bob.cowdery@cgi-europe.com> *** Confidentiality Notice *** Proprietary/Confidential Information belonging to CGI Group Inc. and its affiliates may be contained in this message. If you are not a recipient indicated or intended in this message (or responsible for delivery of this message to such person), or you think for any reason that this message may have been addressed to you in error, you may not use or copy or deliver this message to anyone else. In such case, you should destroy this message and are asked to notify the sender by reply email.

3 2

data smoothing: interpolate.splrep ignores s parameter
by george young Feb. 25, 2004

Feb. 25, 2004

[SciPy-0.2.0_alpha_200.4161, Numeric-23.1, Python 2.3.3, x86 Linux] My goal really is to smooth some noisy measurement data without messing up it's *shape*. My first attempt was 1d splines. I did: from Numeric import * import scipy f = file('rough') data = array([(float(x),float(y)) for (x,y) in [l.split() for l in f]]) rep = scipy.interpolate.splrep(data[:,0], data[:,1], s=s) smooth_y = scipy.interpolate.splev(data[:,0], rep) "s" is supposed to vary the amount of smoothing. For s=0, i get my original data back, as expected. But for all other values, including the recommended range of (m-sqrt(2.0*m)) to (m+sqrt(2.0*m), I get a single, too much smoothed, result. It seems like splrep is not using the "s" value to adjust the splines, except to sense that it's not zero. Splines may not be the right method anyway, since they tend to warp the shape of the curve, and I need to get the data's derivatives. Is there some way to fashion a low pass filter? It seems like fft should be useful here, but I have very little experience with fft's. -- George Young -- "Are the gods not just?" "Oh no, child. What would become of us if they were?" (CSL)

1 0