Thread: Re: cosine taper

Started: 2016-09-26 17:35:01
Last activity: 2016-09-26 17:35:01
Topics: SAC Help
George Helffrich
2016-09-26 17:35:01
Dear All -

I maintain that the COSINE taper is calculated correctly (specifically, by using the sin() function). For justification of this assertion, please refer to the paper,

Harris, F. J. (1978). On the use of windows for harmonic analysis with the discrete Fourier transform. Proc. IEEE 66(1), 51-83.

It deals explicily with the cosine taper, and presents the time-domain coefficients w(n) (n.b. not W(n)) for it. They are calculated through the sin() function, the reason for sin() being the pi/2 phase shift required by reindexing the phase angle from zero at the start of the trace. In particular, see sections IV (eqns. 3a-3c, 4a-4d) and section V(B), equations 25a & 25b.

It seems that the only problem with SAC is the documentation that gives the formula for the time domain taper coefficients, which is perhaps Peter’s point that got lost in the comments on the form of the code. SAC’s code does, however, correctly implement a cosine taper.

If one desires smoother taper behavior at the ends of the window, use the Hann or cos()**n (n>1) tapers (and suffer the spectral leakage around the central lobe of the taper — see Figs. 16-19 of the cited paper).

On 23 Sep 2016, at 22:22, Peter Schmidt <peter.schmidt<at>geo.uu.se> wrote:

Hi George and others

To clarify, the algorithm you posted is exactly the same as used by the C-source I've got. (I merely represented the algorithm by abs(sin(...)) where the abs( ) is used to point out that the taper is based on an even function, which the sin() function is not)

My point here is merely that the manual is misleading in stating that a 1-cos(...) taper is applied when its not.


regP


On 23/09/16 14:57, George Helffrich (ELSI) wrote:
Dear All -

This is not true of the original Fortran code (10.6d or 10.6f). In that source, the taper is calculated from

sin(j*pi/(2*npts)), j = 0, npts

where

npts = w*(e - b)/delta

with e, b and delta having their SAC file header definitions and w is the fractional taper length. Hence, if this is how the C SAC source code behaves, it should be viewed as a bug or as an incompatible change with respect to historical behavior.

On 23 Sep 2016, at 18:40, Peter Schmidt <peter.schmidt<at>geo.uu.se> wrote:

Hi All

Just a small note that perhaps would be reasonable to edit in the SAC manual

According to manual (V101.6a, November 21, 2013) "TAPER TYPE COSINE
WIDTH v" will apply a symmetric taper on format: 1-COS(...) to each end
of the data. However after getting some inconsistent results I had a
look in the source code and found out that the applied taper is really:
ABS(SIN(...)). Perhaps a future edit of the manual could properly state
this.


regP

--
********************************************************************************
Peter Schmidt Tel: +46-18-4712259
Swedish National Seismic Network (SNSN) Mobile: +46-73-3190975
Dept. of Earth Sciences:geophysics e-mail: peter.schmidt<at>geo.uu.se
Uppsala University
Villavagen 16
SE-75236 Uppsala
********************************************************************************

----------------------
SAC Help (http://ds.iris.edu/message-center/topic/sac-help/)

Sent from the IRIS Message Center (http://ds.iris.edu/message-center/)
Update subscription preferences at http://ds.iris.edu/account/profile/


George Helffrich
george<at>elsi.jp


--
********************************************************************************
Peter Schmidt Tel: +46-18-4712259
Swedish National Seismic Network (SNSN) Mobile: +46-73-3190975
Dept. of Earth Sciences:geophysics e-mail: peter.schmidt<at>geo.uu.se
Uppsala University
Villavagen 16
SE-75236 Uppsala
********************************************************************************




George Helffrich
george<at>elsi.jp


  • Peter Schmidt
    2016-09-26 17:14:51
    Hi George

    Thanks for the reference.

    Yes, my point were all about the documentation, but managed to obfuscate
    this.


    regP


    On 26/09/16 03:36, George Helffrich wrote:
    Dear All -

    I maintain that the COSINE taper is calculated correctly (specifically, by using the sin() function). For justification of this assertion, please refer to the paper,

    Harris, F. J. (1978). On the use of windows for harmonic analysis with the discrete Fourier transform. Proc. IEEE 66(1), 51-83.

    It deals explicily with the cosine taper, and presents the time-domain coefficients w(n) (n.b. not W(n)) for it. They are calculated through the sin() function, the reason for sin() being the pi/2 phase shift required by reindexing the phase angle from zero at the start of the trace. In particular, see sections IV (eqns. 3a-3c, 4a-4d) and section V(B), equations 25a & 25b.

    It seems that the only problem with SAC is the documentation that gives the formula for the time domain taper coefficients, which is perhaps Peter’s point that got lost in the comments on the form of the code. SAC’s code does, however, correctly implement a cosine taper.

    If one desires smoother taper behavior at the ends of the window, use the Hann or cos()**n (n>1) tapers (and suffer the spectral leakage around the central lobe of the taper — see Figs. 16-19 of the cited paper).

    On 23 Sep 2016, at 22:22, Peter Schmidt <peter.schmidt<at>geo.uu.se> wrote:

    Hi George and others

    To clarify, the algorithm you posted is exactly the same as used by the C-source I've got. (I merely represented the algorithm by abs(sin(...)) where the abs( ) is used to point out that the taper is based on an even function, which the sin() function is not)

    My point here is merely that the manual is misleading in stating that a 1-cos(...) taper is applied when its not.


    regP


    On 23/09/16 14:57, George Helffrich (ELSI) wrote:
    Dear All -

    This is not true of the original Fortran code (10.6d or 10.6f). In that source, the taper is calculated from

    sin(j*pi/(2*npts)), j = 0, npts

    where

    npts = w*(e - b)/delta

    with e, b and delta having their SAC file header definitions and w is the fractional taper length. Hence, if this is how the C SAC source code behaves, it should be viewed as a bug or as an incompatible change with respect to historical behavior.

    On 23 Sep 2016, at 18:40, Peter Schmidt <peter.schmidt<at>geo.uu.se> wrote:

    Hi All

    Just a small note that perhaps would be reasonable to edit in the SAC manual

    According to manual (V101.6a, November 21, 2013) "TAPER TYPE COSINE
    WIDTH v" will apply a symmetric taper on format: 1-COS(...) to each end
    of the data. However after getting some inconsistent results I had a
    look in the source code and found out that the applied taper is really:
    ABS(SIN(...)). Perhaps a future edit of the manual could properly state
    this.


    regP

    --
    ********************************************************************************
    Peter Schmidt Tel: +46-18-4712259
    Swedish National Seismic Network (SNSN) Mobile: +46-73-3190975
    Dept. of Earth Sciences:geophysics e-mail: peter.schmidt<at>geo.uu.se
    Uppsala University
    Villavagen 16
    SE-75236 Uppsala
    ********************************************************************************

    ----------------------
    SAC Help (http://ds.iris.edu/message-center/topic/sac-help/)

    Sent from the IRIS Message Center (http://ds.iris.edu/message-center/)
    Update subscription preferences at http://ds.iris.edu/account/profile/

    George Helffrich
    george<at>elsi.jp

    --
    ********************************************************************************
    Peter Schmidt Tel: +46-18-4712259
    Swedish National Seismic Network (SNSN) Mobile: +46-73-3190975
    Dept. of Earth Sciences:geophysics e-mail: peter.schmidt<at>geo.uu.se
    Uppsala University
    Villavagen 16
    SE-75236 Uppsala
    ********************************************************************************



    George Helffrich
    george<at>elsi.jp


    ----------------------
    SAC Help (http://ds.iris.edu/message-center/topic/sac-help/)

    Sent from the IRIS Message Center (http://ds.iris.edu/message-center/)
    Update subscription preferences at http://ds.iris.edu/account/profile/

    --
    ********************************************************************************
    Peter Schmidt Tel: +46-18-4712259
    Swedish National Seismic Network (SNSN) Mobile: +46-73-3190975
    Dept. of Earth Sciences:geophysics e-mail: peter.schmidt<at>geo.uu.se
    Uppsala University
    Villavagen 16
    SE-75236 Uppsala
    ********************************************************************************

00:47:08 v.22510d55