graetz23 / JWave

A Discrete Fourier Transform (DFT), a Fast Wavelet Transform (FWT), and a Wavelet Packet Transform (WPT) algorithm in 1-D, 2-D, and 3-D using normalized orthogonal (orthonormal) Haar, Coiflet, Daubechie, Legendre and normalized biorthognal wavelets in Java.
http://en.wikipedia.org/wiki/Wavelet
Other
215 stars 73 forks source link

Flipped signs in all high pass filter coefficients, error in _buildOrthonormalSpace #19

Closed pranathivemuri closed 6 years ago

pranathivemuri commented 6 years ago

Hi!

I found an issue with jwave's wavelet decomposed coefficients having sign differences in detail coefficients when compared with pywavelet coefficients.

So, I digged through your code, and figured that you are using pywavelets mother wavelet decomposition function as low pass filter coefficients from this page (http://wavelets.pybytes.com/wavelet/db10/) and building the rest using _buildOrthonormalSpace here (https://github.com/cscheiblich/JWave/blob/master/src/jwave/transforms/wavelets/Wavelet.java#L104)

The high pass filter coefficients you have build using this function do not turn out to be the same as the high pass filter coefficients for the mother wavelet scaling function as given in pywavelets documentation (http://wavelets.pybytes.com/wavelet/db10/)

Here is a demonstration of error in coefficients output from Daubechies10 as compared to python

Python

In [16]: arr = np.zeros(4)

In [17]: arr[1] = 1

In [18]: coeffs = [0] * 3 ...: for i in range(3): ...: arr, coeffs[i] = pywt.dwt(arr, "db10", "periodization") ...: print(arr, coeffs[i]) ...:
...:
[ 0.30328632 0.40382046] [-0.14391341 0.85102019] [ 0.5] [ 0.07108838] [ 0.70710678] [ -5.07923307e-17]

// Java

/**

level 0 [0.3032863181493727, 0.4038204630371676, 0.14391341375270209, -0.8510201949392565] level 1 [0.49999999999999406, -0.07108837559095565, 0.14391341375270209, -0.8510201949392565]

Error in high pass filter coefficient

In [51]: _scalingDeCom = [-1.3264203002354869e-005, 9.3588670001089845e-005, -0.0001164668549943862, -0.00068585669500468248, 0.0019924052949908499, 0.00139535174699407 ...: 98, -0.010733175482979604, 0.0036065535669883944, 0.033212674058933238, -0.029457536821945671, -0.071394147165860775, 0.093057364603806592, 0.12736934033574265 ...: , -0.19594627437659665, -0.24984642432648865, 0.28117234366042648, 0.68845903945259213, 0.52720118893091983, 0.18817680007762133, 0.026670057900950818] ...: _motherWavelength = len(_scalingDeCom) ...: _waveletDeCom = [0] * _motherWavelength; ...: for i in range(0, _motherWavelength): ...: if( i % 2 == 0 ): ...: _waveletDeCom[ i ] = _scalingDeCom[ ( _motherWavelength - 1 ) - i ]; ...: else: ...: _waveletDeCom[ i ] = -_scalingDeCom[ ( _motherWavelength - 1 ) - i ]; ...:

In [52]: _waveletDeCom Out[52]: [0.026670057900950818, -0.18817680007762133, 0.5272011889309198, -0.6884590394525921, 0.2811723436604265, 0.24984642432648865, -0.19594627437659665, -0.12736934033574265, 0.09305736460380659, 0.07139414716586077, -0.02945753682194567, -0.03321267405893324, 0.0036065535669883944, 0.010733175482979604, 0.0013953517469940798, -0.00199240529499085, -0.0006858566950046825, 0.0001164668549943862, 9.358867000108985e-05, 1.326420300235487e-05]

expected = [-0.026670057900950818 0.18817680007762133 -0.5272011889309198 0.6884590394525921 -0.2811723436604265 -0.24984642432648865 0.19594627437659665 0.12736934033574265 -0.09305736460380659 -0.07139414716586077 0.02945753682194567 0.03321267405893324 -0.0036065535669883944 -0.010733175482979604 -0.0013953517469940798 0.00199240529499085 0.0006858566950046825 -0.0001164668549943862 -9.358867000108985e-05 -1.326420300235487e-05]

Hence the output detail coefficients are differing in signal, since you have used the same flipped high pass filter for building the rest of the coefficients for reconstruction, your inverse transform did not make any difference, but this is an error.

Fix - basically line 112 and 114 signs should be flipped in function _buildOrthonormalSpace, since the first coefficient in high pass filter should have a negative sign -(https://github.com/cscheiblich/JWave/blob/master/src/jwave/transforms/wavelets/Wavelet.java#L104)

protected void _buildOrthonormalSpace( ) { // building wavelet as orthogonal (orthonormal) space from // scaling coefficients (low pass filter). Have a look into // Alfred Haar's wavelet or the Daubechies Wavelet with 2 // vanishing moments for understanding what is done here. ;-) _waveletDeCom = new double[ _motherWavelength ]; for( int i = 0; i < _motherWavelength; i++ ) if( i % 2 == 0 ) _waveletDeCom[ i ] = -_scalingDeCom[ ( _motherWavelength - 1 ) - i ]; else _waveletDeCom[ i ] = _scalingDeCom[ ( _motherWavelength - 1 ) - i ]; // Copy to reconstruction filters due to orthogonality (orthonormality)! _scalingReCon = new double[ _motherWavelength ]; _waveletReCon = new double[ _motherWavelength ]; for( int i = 0; i < _motherWavelength; i++ ) { _scalingReCon[ i ] = _scalingDeCom[ i ]; _waveletReCon[ i ] = _waveletDeCom[ i ]; } // i } // _buildOrthonormalSpace

cscheiblich commented 6 years ago

Hi there,

acutally there is no error! :-)

Some of the base vector pairs of the multidimensional orthogonal (orthonormal) space (wavelet or scaling) are differently ordered in pyWavelets (flipping signs; see your comment). From filtering effect and /or quality JWave makes no difference! The resulting coeffiecients in Hilbert space are just ordered and sampled differently (entropy, etc. are the same) - that's all.

By the way, I reorded those to harmonize the buildSpace method with the other wavelets of JWave. ;-)

BR Chris