Hi Tim, Brilliant! Many thanks... I think this is exactly what I need, I owe you a beer (or other beverage of your choice).
I'm now going to lock myself in the basement until I can work out an implementation of this for my use-case :) /Carl On Tue, Sep 3, 2013 at 9:05 PM, Cera, Tim <t...@cerazone.net> wrote: > > I am trying to take the rfft of a numpy array, like this: > >>>> my_rfft = numpy.fft.rfft(my_numpy_array) > > > > and replace the amplitudes that can be obtained with: > > > >>>> my_amplitudes = numpy.abs(my_rfft) > > > > with amplitudes from an arbitrary numpy array's rFFT, which is to then be > > converted back using numpy.fft.irfft . Alternately, some future plans > will > > involve having to modify individual array element amplitudes directly > based > > on other parameters. I would think that modifying and re-synthesizing > > signals using FFT is a fairly common use-case, but my attempts at > Googling > > example code have been fruitless. > > I have FFT transform filter in my tidal analysis package. See > > http://sourceforge.net/apps/mediawiki/tappy/index.php?title=CompareTidalFilters > for a comparison and short description. > > See my function below. My earlier self made some poor variable name > choices. The 'low_bound' variable is actually where frequencies > greater are set to zero ('factor[freq > low_bound] = 0.0'), then > factor is ramped from 0 at 'low_bound' to 1 at 'high_bound'. To > filter out tidal signals if your water elevations are hourly then > 'low_bound' = 1/30.0 and 'high_bound' = 1/40.0. Having this gradual > change in the frequency domain rather than an abrupt change makes a > better filter. > > def fft_lowpass(nelevation, low_bound, high_bound): > """ Performs a low pass filter on the nelevation series. > low_bound and high_bound specifies the boundary of the filter. > """ > import numpy.fft as F > if len(nelevation) % 2: > result = F.rfft(nelevation, len(nelevation)) > else: > result = F.rfft(nelevation) > freq = F.fftfreq(len(nelevation))[:len(nelevation)/2] > factor = np.ones_like(result) > factor[freq > low_bound] = 0.0 > > sl = np.logical_and(high_bound < freq, freq < low_bound) > > a = factor[sl] > # Create float array of required length and reverse > a = np.arange(len(a) + 2).astype(float)[::-1] > > # Ramp from 1 to 0 exclusive > a = (a/a[0])[1:-1] > > # Insert ramp into factor > factor[sl] = a > > result = result * factor > print 'result=', len(result) > relevation = F.irfft(result, len(nelevation)) > print 'result=', len(relevation) > return relevation > > > Kindest regards, > Tim > _______________________________________________ > NumPy-Discussion mailing list > NumPy-Discussion@scipy.org > http://mail.scipy.org/mailman/listinfo/numpy-discussion >
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
