HOW-TO-Decimate

The Art of Data Decimation

If you have to down-sample 1-s data to e.g. 10-min, the required low-pass filter will be excessively long and convolution costly.
Doing it by Fourier transform does not work if data points are missing. A possible solution is

Down-sampling in a cascade. An example:

TSF EDIT 600REPAIR L<=10FILTER D:WD KAIS 40 2.1 0.0 0.05 0.5 0.0DECIMATE 10 ${MOVE} 0FILTER D:WD KAIS 58 2.1 0 0.0333 1.0 0.0DECIMATE 60 ${MOVEKWD30} 27END

One can show that the 10+60 sequence is optimum for 1-to-600 decimation
What is the harmonic response of the cascade? We can easily compute it for the second stage. The objective is to preserve (and precisely determine) the harmonic response to tides.

So let's make a test, 3 days of 1-s data and an M2 signal. This is a segment of prf.tse:

TSF EDIT SINSIN 1.0 ${SINFRQ:[0.0805114]} 0.0FILTER D:WD KAIS 40 2.1 0 0.05 1.0 0.0DECIMATE 10 ${MOVE} 0FILTER D:WD KAIS 58 2.1 0 0.0333 1.0 0.0FILTHR TAPRL=31 R t/urtap.prl]DECIMATE 30 ${MOVEKWD30} 27END

tslist _259200 -B -N -rs1. -Eprf.tse,SIN -qm -I

The -qm option gives us the min and max of the output, theFILTHR command the response of the second stage to a number of tides.
( MOVE and MOVEKWD were unset.)

Result:
<TsfEdi>>> M2 1.9322740 [cyc/d] 12.4206 [h] :: 0.999726 0.000000 <Main-->>> Min = -9.9993E-01 at 279 D T = 2013 10 25 22 34 50 000 <Main-->>> Max = 9.9985E-01 at 18 D T = 2013 10 24 03 04 50 000
You may give some grace since the output is sampled at 600s. If the last DECIMATE command is avoided, we get instead
<Main-->>> Min = -9.9993E-01 at 7764 D T = 2013 10 24 21 44 10 000 <Main-->>> Max = 9.9993E-01 at 23414 D T = 2013 10 26 17 12 30 000
Nice! So you don't need the response of the first stage.
Next example is for an additional sine that will be aliased to M2.
The sampling rate is 0.05 Hz. A slightly lower frequency will be aliased more strongly than a slightly higher, so

`setenv FRQ = `calc -f "%15.10f" '3600*0.05-0.0805114'

and synthesize the following signal in prf.tse,SIN :

SIN 1.0 0.0805114 0.0
SIN 1.0 ${FRQ:[179.9194886000]} 0.0

Result:
<Main-->>> Min = -9.9993E-01 at 7764 D T = 2013 10 24 21 44 10 000 <Main-->>> Max = 9.9993E-01 at 14471 D T = 2013 10 25 16 22 00 000
No real difference (except that the date of the occurrence of the maximum has changed).

The spectrum of the filter looks convincing:

sasm02 -f WD KAIS 40 2.1 0 0.05 0.5 0.0

and with the ?f 0.0499 request we find -54 dB at a frequency slightly below the one aliased to M2.

The objective should be that the attenuation at frequencies aliasing to M3-like frequencies is below 80 dB. Here's a better filter
(greater design-parameter, slightly lower passband corner, same length)

sasm02 -f WD KAIS 40 2.5 0 0.035 0.5 0.0

A slightly narrower low-pass is here:
sasm02 -d8192 -f WD KAIS 40 2.5 0 0.01 0.5 0.0
We should hone the edge so that the second node occurs exactly at 0.05 Hz. A good guess is 0.013.
<TsfEdi>>> M2 1.9322740 [cyc/h] 12.4206 [h] :: 0.999931 0.000000 <Main-->>> Min = -9.999302E-01 at 7764 D T = 2013 10 24 21 44 10 000 <Main-->>> Max = 9.999302E-01 at 23414 D T = 2013 10 26 17 12 30 000This is good to 1:10⁵, enough.

Summary: We have a filter that passes M2 and stops a signal that would be aliased to M2 at -100 dB.

We should not try to adjust the filter with sasm02. tslist/tsfedit FILTHR is the preferred method.

So far we looked at stage one of the cascade. Now stage 2.
Here, we do not need a flat passband, since the transfer function can be used in the tide analysis procedure (urtapt) to compensate for the attenuation.
If we require ~ -80 dB for the stopband, the length of the filter will have to be > 240.

sasm02 -d8192 -f WD KAIS 240 2.5 0.0 0.0025 0.5 0.0

The Nyquist of the 600 s resampling is at 0.0083333 Hz (green number, reads -80.765 dB).
The first point of -80 dB is at
0.00775 Hz, 0.93 of the Nyquist.
The damping at M2 is
<TsfEdi>>> M2 1.9322740 [cyc/h] 12.4206 [h] :: 0.996095 0.000000

TSF EDIT from 1s to 600s with a two-stage cascade

TSF EDIT C600S; first record out will be at 00:50:00 FILTER D:WD KAIS 40 2.5 0 ${UBE:[0.013]} 0.5 0.0 DECIMATE 10 ${MOVE:[560]} 0 FILTER D:WD KAIS 240 2.5 0.0 0.0025 0.5 0.0 DECIMATE 60 ${MOVEKWD60:[0]} END

TSF EDIT from 1s to 1h with a three-stage cascade

TSF EDIT C1H; first record out will be at 05:00:00FILTER D:WD KAIS 48 2.5 0.0 0.039 0.5 0.0DECIMATE 5 3272FILTER D:WD KAIS 56 2.5 0 0.008333333 0.1 0.0DECIMATE 12 0FILTER D:WD KAIS 240 2.5 0 0.00833333 0.5 0.0DECIMATE 60 0END

The filter series of the last stage have been saved in ~/TD/f (linked ~/Ttide/SCG/f ):

-rw-rw-r-- 1 hgs hgs 2033 Oct 31 00:35 f/WD-KAIS-240-2.5-0.0025-0.5-600s.fs-rw-rw-r-- 1 hgs hgs 2033 Oct 31 00:39 f/WD-KAIS-240-2.5-0.00833-0.5-1h.fs

There are plots showing pulse responses of the decimation schemes. If the pulse length reads 0, the pulse is a single spike. In all other cases, the pulse has a zero-mean. The lengths shown mean full pulse length (the parameter for generation a pulse designate the number of positive off-centre ordinates). The pulse generator is here: ~/sas/p/pulsegens.f , double precision function pulsegen (t,scale,lhwidth).

do i=-lhwidth,lhwidth xt=dfloat(i)/hw xv=(dcos(dpi/2.d0*xt))**2 & *( (dcos(dpi*xt))**2 & - (cos(dpi*(dabs(xt)-0.5d0)))**2 ) x(i)=xt enddo


	Obs! that the longest pulse extrudes beyond the decimation filter length of the last stage.

MATHS
A short and steep filter will have to have a length of 8 D +1 where D is the decimation factor. For a 1:10 decimation we need 81 coefficients (it's a symmetric filter with 40 off-centre coefficients). We will call this factor 8 the decimation filter scale.

The number of operations is proportional to D. In a cascade the number of operations is

N_op = N (8 D₁ +1 ) + N/D₁ (8 D₂ + 1 )
D₁ D₂ = D_tot

Thus,

N_op = N (8 D₁ + 1 ) + N/D₁ (8 D_tot/D₁ + 1 )

In the number of operations (per sample of the input series times the filter scale)

N_op / 8N = D₁ + 1/8 + D_tot/D₁² + 1/(8 D₁)

we can neglect the last term. The minimum is found at

0 = 1 - 2 D_tot/D₁³

D₁ = cubicroot( 2 D_tot )

Examples:

D_tot = 600 => D₁ = 10.63

D_tot = 3600 => D₁ = 19.31

.bye