Re: Question on the low pass filter?

From: Dennis Shea <shea_at_nyahnyahspammersnyahnyah>
Date: Tue, 26 Dec 2006 11:43:27 -0700 (MST)

>Dear NCL user:

>Hi, I want to do a low pass filter (filtering signal
> with period less than 10a) for a 100a time series.
> After see some tips on the forum , I still have some qestions.
> Using "ezfftf" and "ezfftb" :
> x = random_uniform(-1,1,100) ; using a random series as example
> fc = ezfftf (x)
> fc (:,5:49) = 0.0 ; Is it right for filtering the
> ; period longer than 10a?
> xband = ezzftb(fc) ; the result
> Using filwgts_lancos :
> I want to know how to define the nwt and fca
> to get a similar low pass filter as listed above ?

I am not sure what the "a" in "10a" represents.

Filtering via FFT, simplistic explanation:

  [1] Truncating a fourier series as you have done above
      results in Gibbs phenomena. When the series is long,
      there will be a 9% over estimate at the truncated
      frequency [fc] and a 5% under estimate at 2*fc.
      If you look at some text book or google on Gibbs
      phenomena you will note that there will *always*
      be some "ripples"

  [2] FFTs should be performed on periodic sequences and
      time series are not periodic.
      [a] Remove the original series mean. [xMean]
          The mean can be added back in later.
      [b] Taper the series via

          This results is some loss of variance but will
          minimize 'spectral leakage'.
  [3] You must weight the FFT coefficients in a transition
                                   ideal low pass
             1.0 i---------
                  i |
                  i |
                  i | This results in Gibbs
                  i |
             0.0 i________|______________________________
                 0.0 fca 0.5

                                   transition zone low pass
             1.0 i---------
                  i \
                  i \ the \ indicate the transition
                  i \
                  i \
             0.0 i_____________\_________________________
                 0.0 fca 0.5

      For illustration
                 fc(:,5) = 0.75*fc(:,5)
                 fc(:,6) = 0.50*fc(:,6)
                 fc(:,7) = 0.25*fc(:,7)
                 fc(:,8:)= 0.0
      large transition zone: small ripples
      small " " : large ripples but less than Gibbs
  [4] Use ezfftb to go back to the time domain.
  [5] If you have many series or *very* long
      series then the the FFT has a speed advantage
      over classic weighting of values in the time
  [6] Some people mistakingly think that since the FFT
      returns the original number of points [N, here 100]
      that the FFT is better than applying weights
      to the series in the time domain. After all, in addition
      to being faster, the latter results in a few points
      at the beginning and end being set to missing.
      Not so. There is no magic. At least with time domain
      weighting, the user can be assured that the resulting
      series has consistent frequency response across all values.
      I would be skeptical of the first/last few points
      generated via ezfftb in subsequent analysis. How many
      points? I am not aware of any 'formula' that give the
      answer. One approach would be to use lancos weighting
      to a series; then, use trial and error to vary the
      transition zone weighting to the fft coefficients.
      Plot the test series.
> Using filwgts_lancos :
> I want to know how to define the nwt and fca
> to get a similar low pass filter as listed above ?

df = 1./100
fca = 10*df ; for '10a'

nwt = ; this is art
              ; there is a trade off between the ideal
              ; and what is practically possible.
              ; The smaller nwt the fewer points lost BUT
              ; the poorer response of the filter; the larger
              ; nwts, the more points lost but a better filter.
              ; Again, there is no magic. The fft has the
              ; same issues but in a different form.

I have a script that I will send you off line.
There will be *no support* for this offline script.good luck

good luck

ncl-talk mailing list
Received on Tue Dec 26 2006 - 11:43:27 MST

This archive was generated by hypermail 2.2.0 : Thu Dec 28 2006 - 09:40:18 MST