<jupyter_start><jupyter_text>Check whether the Silver and Chan (1991) coefficients for their implementation of Parseval's theorem is correct for digitised data.<jupyter_code>%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import math
from scipy import signal
from scipy import stats<jupyter_output><empty_output><jupyter_text>Parsevals theorem when applied to discrete Fourier Transform looks like this.$\sum _{n=0}^{N-1}|x[n]|^{2}={\frac {1}{N}}\sum _{k=0}^{N-1}|X[k]|^{2}$Source: https://en.wikipedia.org/wiki/Parseval%27s_theorem<jupyter_code># check Parseval's theorem holds numerically 
nsamps=1000

# window
w = signal.tukey(nsamps,0.1)

a = np.random.normal(0,1,nsamps) * w
A = np.fft.fft(a)
b = (1/np.sqrt(2*np.pi))*(signal.gaussian(nsamps,10))
B = np.fft.fft(b)
c = np.convolve(a,b,'same')
C = np.fft.fft(c)

# signal c is convolution of Gaussian noise (a) with a Gaussian wavelet (b)
# C is the fourier transform of c.

sumt = np.sum(c**2)
sumf = np.sum(np.abs(C)**2)/nsamps

print('time domain',sumt)
print('fourier domain',sumf)
print('difference',np.abs(sumt-sumf))
print('percent', (np.abs(sumt-sumf)/sumt)*100)<jupyter_output>time domain 2778.93808225
fourier domain 2778.93808225
difference 9.09494701773e-13
percent 3.27281384059e-14<jupyter_text>Furthermore by the convolution theorem: C = A \* B.And therefore sum(C^2) = sum(A^2 * B^2)<jupyter_code>AB = A * B
ab = np.fft.ifft(AB)
plt.plot(np.roll(ab,500))
plt.plot(c)

sumAB = np.sum(np.abs(A**2*B**2))/nsamps
print('sum A*B',sumAB)
print('difference',np.abs(sumt-sumAB))
print('percent',(np.abs(sumt-sumAB)/sumt)*100)<jupyter_output>/Users/glyjw/anaconda/envs/py35/lib/python3.5/site-packages/numpy/core/numeric.py:531: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)<jupyter_text>Parsevals theorem as applied in Silver and Chan (and Walsh).$\sum _{n=0}^{N-1}|x[n]|^{2}={\frac {1}{N}}\sum _{k=1}^{N-2}|X[k]|^{2}+\frac{1}_{2}\sum|X[0,N-1]|$Source: https://en.wikipedia.org/wiki/Parseval%27s_theorem<jupyter_code>def ndf(y,taper=True,detrend=True):
    """
    Uses the improvement found by Walsh et al (2013).
    By default will detrend data to ensure zero mean
    and will taper edges using a Tukey filter affecting amplitudes of 5% of data at edges
    """

    if taper is True:
        y = y * signal.tukey(y.size,0.05)
        
    if detrend is True:
        # ensure no trend on the noise trace
        y = signal.detrend(y)

  
    Y = np.fft.fft(y)
    amp = np.absolute(Y)
    
    # estimate E2 and E4 following Walsh et al (2013)
    a = np.ones(Y.size)
    a[0] = a[-1] = 0.5
    E2 = np.sum( a * amp**2)
    E4 = (np.sum( (4 * a**2 / 3) * amp**4))
    
    ndf = 2 * ( 2 * E2**2 / E4 - 1 )
    
    return ndf
    
def ndf2(y,taper=True,detrend=True):
    """
    
    """

    if taper is True:
        y = y * signal.tukey(y.size,0.05)
        
    if detrend is True:
        # ensure no trend on the noise trace
        y = signal.detrend(y)

  
    Y = np.fft.fft(y)
    amp = np.absolute(Y)**2
    
    E2 = np.sum(amp**2)
    E4 = (np.sum( (4/3) * amp**4))
    
    ndf = 2 * ( 2 * E2**2 / E4 - 1 )
    
    return ndf
print(ndf(c))
print(ndf2(c))
stats.moment(c,moment=4)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16150576806314512, 'fraction_numerical': 0.04948391013964785, 'mean_word_length': 3.6343178621659633, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 15, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '230159', 'license': 'MIT', 'n_tokens_mistral': 1415, 'n_tokens_neox': 1301, 'n_words': 299, 'path': 'devel/Parseval.ipynb', 'repo_name': 'kaushalkapadiya/splitwavepy', 'size': '3294'}