From: Dennis Shea <shea_at_nyahnyahspammersnyahnyah>

Date: Tue, 26 Dec 2006 11:43:27 -0700 (MST)

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

http://www.ncl.ucar.edu/Document/Functions/Built-in/taper.shtml

This results is some loss of variance but will

minimize 'spectral leakage'.

[3] You must weight the FFT coefficients in a transition

zone.

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

domain.

[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.

See:

http://www.ncl.ucar.edu/Document/Functions/Built-in/filwgts_lancos.shtml

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

ncl-talk_at_ucar.edu

http://mailman.ucar.edu/mailman/listinfo/ncl-talk

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
*