001: package JSci.maths.wavelet.daubechies4;
002:
003: import JSci.maths.wavelet.*;
004: import JSci.maths.*;
005:
006: /******************************************
007: * Daubechies wavelets adapted to the
008: * interval by Meyer. Thanks to Pierre Vial
009: * for the filters.
010: * @author Daniel Lemire
011: *****************************************/
012: public final class Daubechies4 extends Multiresolution implements
013: Filter, NumericalConstants {
014: protected final static int filtretype = 6;
015: protected final static int minlength = 12;
016:
017: /****************************************
018: * This method is used to compute
019: * how the number of scaling functions
020: * changes from on scale to the other.
021: * Basically, if you have k scaling
022: * function and a Filter of type t, you'll
023: * have 2*k+t scaling functions at the
024: * next scale (dyadic case).
025: * Notice that this method assumes
026: * that one is working with the dyadic
027: * grid while the method "previousDimension"
028: * define in the interface "Filter" doesn't.
029: ******************************************/
030: public int getFilterType() {
031: return (filtretype);
032: }
033:
034: public MultiscaleFunction primaryScaling(int n0, int k) {
035: return (new Scaling4(n0, k));
036: }
037:
038: public MultiscaleFunction dualScaling(int n0, int k) {
039: return (new Scaling4(n0, k));
040: }
041:
042: public MultiscaleFunction primaryWavelet(int n0, int k) {
043: return (new Wavelet4(n0, k));
044: }
045:
046: public MultiscaleFunction dualWavelet(int n0, int k) {
047: return (new Wavelet4(n0, k));
048: }
049:
050: final static double[] vg = { -0.107148901418, -0.0419109651251,
051: 0.703739068656, 1.13665824341, 0.421234534204,
052: -0.140317624179, -0.0178247014417, 0.045570345896 };
053:
054: final static double[] v0temp = { 0.7983434920E+00, 0.6022023488E+00 };
055: final static double[] v1temp = { -0.3918024327E-01,
056: 0.5194149822E-01, -0.4817281609E+00, 0.8739021503E+00 };
057: final static double[] v2temp = { 0.1774707150E-01,
058: -0.2352740580E-01, -0.1232594861E+00, -0.6575127688E-01,
059: -0.9620570014E-01, 0.9850684416E+00 };
060: final static double[] v3temp = { -0.2636405192E-01,
061: 0.3495099166E-01, 0.8114147375E+00, 0.4440233637E+00,
062: 0.3192581817E+00, 0.1636579832E+00, -0.4282797155E-01,
063: 0.1094933054E+00 };
064: final static double[] v4temp = { -0.1670338745E-01,
065: 0.2214378721E-01, -0.1643714751E-01, -0.1112580065E-01,
066: 0.2995602574E+00, 0.2728668922E-01, 0.8472064764E+00,
067: -0.4270166998E+00, -0.3309408518E-01, 0.8460780753E-01 };
068: final static double[] v5temp = { 0.2727915769E-02,
069: -0.3616415322E-02, -0.5206157868E-01, -0.2836107693E-01,
070: -0.4413123462E-01, -0.1285294872E-01, 0.4543141690E+00,
071: 0.8282235028E+00, 0.3000539798E+00, -0.1037443976E+00,
072: -0.1262470890E-01, 0.3227612835E-01 };
073:
074: final static double[] vd0temp = { 0.7629809303E+00,
075: -0.6464209928E+00 };
076: final static double[] vd1temp = { 0.1555526564E+00,
077: 0.1836012627E+00, 0.4620817399E+00, -0.8535657052E+00 };
078: final static double[] vd2temp = { 0.3793246643E+00,
079: 0.4477229057E+00, 0.4284467089E+00, 0.3973740378E+00,
080: -0.2021221018E+00, -0.5228106220E+00 };
081: final static double[] vd3temp = { 0.2385999808E+00,
082: 0.2816233343E+00, 0.1056438723E+00, 0.1612498770E+00,
083: 0.8548427132E+00, 0.2929411663E+00, -0.3647382801E-01,
084: -0.9324840384E-01 };
085: final static double[] vd4temp = { 0.6526723701E-01,
086: 0.7703595299E-01, 0.7744666349E-01, 0.7039069048E-01,
087: -0.6529437593E-01, 0.2555397028E+00, 0.8099281093E+00,
088: 0.4965300820E+00, -0.2995738718E-01, -0.7658857572E-01 };
089: final static double[] vd5temp = { 0.1517778948E-02,
090: 0.1791458518E-02, -0.3127686151E-02, -0.1031248163E-02,
091: 0.3237561439E-01, -0.1322822647E-01, -0.9898430026E-01,
092: 0.2979273659E+00, 0.8037308261E+00, 0.4975940939E+00,
093: -0.2963560969E-01, -0.7576592454E-01 };
094: final static double[] v0 = ArrayMath.scalarMultiply(SQRT2, v0temp);
095: final static double[] v1 = ArrayMath.scalarMultiply(SQRT2, v1temp);
096: final static double[] v2 = ArrayMath.scalarMultiply(SQRT2, v2temp);
097: final static double[] v3 = ArrayMath.scalarMultiply(SQRT2, v3temp);
098: final static double[] v4 = ArrayMath.scalarMultiply(SQRT2, v4temp);
099: final static double[] v5 = ArrayMath.scalarMultiply(SQRT2, v5temp);
100:
101: final static double[] vd0 = ArrayMath.invert(ArrayMath
102: .scalarMultiply(SQRT2, vd0temp));
103: final static double[] vd1 = ArrayMath.invert(ArrayMath
104: .scalarMultiply(SQRT2, vd1temp));
105: final static double[] vd2 = ArrayMath.invert(ArrayMath
106: .scalarMultiply(SQRT2, vd2temp));
107: final static double[] vd3 = ArrayMath.invert(ArrayMath
108: .scalarMultiply(SQRT2, vd3temp));
109: final static double[] vd4 = ArrayMath.invert(ArrayMath
110: .scalarMultiply(SQRT2, vd4temp));
111: final static double[] vd5 = ArrayMath.invert(ArrayMath
112: .scalarMultiply(SQRT2, vd5temp));
113: /********************************************
114: * On définit ici le filtre comme tel par le
115: * vecteur phvg (filtre passe-haut).
116: *********************************************/
117: final static double[] vgtemp = ArrayMath.scalarMultiply(
118: 1.0 / SQRT2, vg);
119: final static double[] phvg = WaveletMath.lowToHigh(vgtemp);
120: final static double[] phv0 = { 0.5979027428E+00, -0.7926434769E+00,
121: -0.1659403671E-01, 0.6477069526E-01, 0.9713044594E-01,
122: -0.1797030610E-01, -0.3192898087E-03, 0.8162911886E-03 };
123: final static double[] phv1 = { 0.4823971249E-01, -0.6395169431E-01,
124: 0.3010034664E+00, 0.1718883936E+00, -0.8873256413E+00,
125: -0.3991915695E-01, 0.2462565991E+00, -0.1524149055E+00,
126: -0.9080945357E-02, 0.2321619929E-01 };
127: final static double[] phv2 = { -0.1162436086E-02, 0.1541048928E-02,
128: 0.2218479707E-01, 0.1208539488E-01, 0.1880547055E-01,
129: 0.5476976811E-02, -0.8114432093E-01, -0.3089428579E+00,
130: 0.8028339176E+00, -0.4957693566E+00, -0.2962669767E-01,
131: 0.7574314021E-01 };
132: final static double[] phvd0temp = { -0.4071236735E+00,
133: -0.4805345164E+00, 0.7323866385E+00, 0.2189246571E+00,
134: 0.1377492261E+00, 0.6432723244E-02, -0.5725128723E-03,
135: -0.1463677232E-02 };
136: final static double[] phvd1temp = { -0.1510687974E+00,
137: -0.1783088929E+00, -0.2220387683E+00, -0.1860863787E+00,
138: 0.4453190824E+00, -0.7578721319E+00, 0.2811170011E+00,
139: 0.9317065817E-01, -0.1190456353E-01, -0.3043501623E-01 };
140: final static double[] phvd2temp = { -0.3568796396E-02,
141: -0.4212306878E-02, 0.7354216552E-02, 0.2424802855E-02,
142: -0.7612569411E-01, 0.3110390153E-01, 0.4970737955E+00,
143: -0.8039165747E+00, 0.2978757587E+00, 0.9927560396E-01,
144: -0.1260377436E-01, -0.3222260742E-01 };
145: final static double[] phvd0 = ArrayMath.invert(phvd0temp);
146: final static double[] phvd1 = ArrayMath.invert(phvd1temp);
147: final static double[] phvd2 = ArrayMath.invert(phvd2temp);
148:
149: /****************************************
150: * This method return the number of "scaling"
151: * functions at the previous scale given a
152: * number of scaling functions. The answer
153: * is always smaller than the provided value
154: * (about half since this is a dyadic
155: * implementation). This relates to the same idea
156: * as the "Filter type". It is used by
157: * the interface "Filter".
158: *****************************************/
159: public int previousDimension(int k) {
160: return (Cascades.previousDimension(filtretype, k));
161:
162: }
163:
164: public Daubechies4() {
165: }
166:
167: /****************************************
168: * This is the implementation of the lowpass
169: * Filter. It is used by the interface
170: * "Filter". Lowpass filters are normalized
171: * so that they preserve constants away from
172: * the boundaries.
173: *****************************************/
174: public double[] lowpass(double[] v, double[] param) {
175: return (lowpass(v));
176: }
177:
178: /****************************************
179: * This is the implementation of the highpass
180: * Filter. It is used by the interface
181: * "Filter". Highpass filters are normalized
182: * in order to get L2 orthonormality of the
183: * resulting wavelets (when it applies).
184: * See the class DiscreteHilbertSpace for
185: * an implementation of the L2 integration.
186: *****************************************/
187: public double[] highpass(double[] v, double[] param) {
188: return (highpass(v));
189: }
190:
191: /****************************************
192: * This is the implementation of the lowpass
193: * Filter. It is used by the interface
194: * "Filter". Lowpass filters are normalized
195: * so that they preserve constants away from
196: * the boundaries.
197: *****************************************/
198: public double[] lowpass(double[] gete) {
199: if (gete.length < minlength) {
200: throw new IllegalScalingException(
201: "The array is not long enough : " + gete.length
202: + " < " + minlength);
203: }
204: double[] sortie = new double[2 * gete.length - filtretype];
205: int dl0 = gete.length - 1;
206: for (int k = 6; k <= dl0 - 6; k++) {
207: for (int L = -4; L < 4; L++) {
208: sortie[2 * k + L - 2] += vg[L + 4] * gete[k];
209: }
210: }
211: sortie = ArrayMath.add(sortie, gete[0], v0, 0);
212: sortie = ArrayMath.add(sortie, gete[1], v1, 0);
213: sortie = ArrayMath.add(sortie, gete[2], v2, 0);
214: sortie = ArrayMath.add(sortie, gete[3], v3, 0);
215: sortie = ArrayMath.add(sortie, gete[4], v4, 0);
216: sortie = ArrayMath.add(sortie, gete[5], v5, 0);
217: int p0 = sortie.length - vd0.length;
218: int p1 = sortie.length - vd1.length;
219: int p2 = sortie.length - vd2.length;
220: int p3 = sortie.length - vd3.length;
221: int p4 = sortie.length - vd4.length;
222: int p5 = sortie.length - vd5.length;
223: sortie = ArrayMath.add(sortie, gete[dl0], vd0, p0);
224: sortie = ArrayMath.add(sortie, gete[dl0 - 1], vd1, p1);
225: sortie = ArrayMath.add(sortie, gete[dl0 - 2], vd2, p2);
226: sortie = ArrayMath.add(sortie, gete[dl0 - 3], vd3, p3);
227: sortie = ArrayMath.add(sortie, gete[dl0 - 4], vd4, p4);
228: sortie = ArrayMath.add(sortie, gete[dl0 - 5], vd5, p5);
229: return (sortie);
230: }
231:
232: /****************************************
233: * This is the implementation of the highpass
234: * Filter. It is used by the interface
235: * "Filter". Highpass filters are normalized
236: * in order to get L2 orthonormality of the
237: * resulting wavelets (when it applies).
238: * See the class DiscreteHilbertSpace for
239: * an implementation of the L2 integration.
240: *****************************************/
241: public double[] highpass(double[] gete) {
242: double[] sortie = new double[2 * gete.length + filtretype];
243: int dl0 = gete.length - 1;
244: for (int k = 3; k <= dl0 - 3; k++) {
245: for (int L = -4; L < 4; L++) {
246: sortie[2 * k + L + 4] += phvg[L + 4] * gete[k];
247: }
248: }
249: sortie = ArrayMath.add(sortie, gete[0], phv0, 0);
250: int p0 = sortie.length - phvd0.length;
251: sortie = ArrayMath.add(sortie, gete[dl0], phvd0, p0);
252: sortie = ArrayMath.add(sortie, gete[1], phv1, 0);
253: int p1 = sortie.length - phvd1.length;
254: sortie = ArrayMath.add(sortie, gete[dl0 - 1], phvd1, p1);
255: sortie = ArrayMath.add(sortie, gete[2], phv2, 0);
256: int p2 = sortie.length - phvd2.length;
257: sortie = ArrayMath.add(sortie, gete[dl0 - 2], phvd2, p2);
258:
259: return (sortie);
260:
261: }
262:
263: public double[] evalScaling(int n0, int k, int j1) {
264: return (Cascades.evalScaling(this , n0, j1, k));
265: }
266:
267: public double[] evalWavelet(int n0, int k, int j1) {
268: return (Cascades.evalWavelet(this, filtretype, n0, j1, k));
269: }
270:
271: }
|