001    package util.wavelet;
002    import java.util.List;
003    import static java.lang.Math.*;
004    
005    /** A wrapper for an array of signal values and boundary conditions check.
006     * @author John Talbot
007     */
008    public class SimpleSignal {
009           
010        /** Signal value array. */
011        float[] values;
012    
013        /** 
014         * Extrapolation methodology for computing signal values for out of range positions. 
015         * Default is mirror boundary.
016         */
017        BoundaryCondition boundaryCondition = BoundaryCondition.Mirror;
018    
019        /** 
020         * Construct a signal using the given array of signal values and default boundary conditions.
021         * @param values an array of signal values
022         */
023        public SimpleSignal(float[] values) {
024            this.values = values;
025        }
026    
027        /** 
028         * Construct a signal using the given array of signal values and boundary conditions.
029         * @param values an array of signal values
030         * @param boundaryCondition the boundary condition
031         */
032        public SimpleSignal(float[] values, BoundaryCondition boundaryCondition) {
033            this(values);
034            this.boundaryCondition = boundaryCondition;
035        }
036    
037        /** 
038         * Copy constructor.
039         * @param signal to be cloned
040         */
041        public SimpleSignal(SimpleSignal signal) {
042            this.boundaryCondition = signal.getBoundaryCondition();
043            values = new float[signal.size()];
044            System.arraycopy(signal.values, 0, values, 0, signal.size());
045        }
046    
047        /** 
048         * Construct a zero filled signal using the size and default boundary conditions.
049         * @param size the size of the signal
050         */
051        public SimpleSignal(int size) {
052            this(new float[size]);
053        }
054    
055        /** 
056         * Construct a zero filled signal using the size and a boundary condition.
057         * @param size the size of the signal
058         * @param boundaryCondition the boundary condition
059         */
060        public SimpleSignal(int size, BoundaryCondition boundaryCondition) {
061            this(size);
062            this.boundaryCondition = boundaryCondition;
063        }
064    
065        /** 
066         * Convolve this signal with a kernel. Not efficient enough for atrous DWT.
067         * For symmetric convolutions use only odd sized kernels.
068         * @param kernel the convolution kernel
069         */
070        public SimpleSignal convolve(SimpleSignal kernel) {
071            int offset = kernel.size() / 2;
072            for (int i = 0; i < size(); i++) {
073                float sum = 0;
074                for (int j = 0; j < kernel.size(); j++)
075                   sum += getValue(i + j - offset) * kernel.getValue(j);
076                values[i] = sum;
077            }        
078            return this;
079        }
080    
081        /** 
082         * Get the boundary condition (extrapolation methodology) for this signal.
083         * @return the boundary condition;
084         */
085        public BoundaryCondition getBoundaryCondition() {
086            return boundaryCondition;
087        }
088    
089        /** 
090         * Set the boundary condition (extrapolation methodology) for this signal.
091         * @param boundaryCondition  the boundary condition
092         */
093        public void setBoundaryCondition(BoundaryCondition boundaryCondition) {
094            this.boundaryCondition = boundaryCondition;
095        }
096           
097        /** 
098         * Get the size of this signal.
099         * @return size of the signal
100         */
101        public int size() {
102            return values.length;
103        }
104    
105        /** 
106         * Get the signal value at the given position and deal with boundary conditions when position is out of range.
107         * @param    position    the position
108         * @return               value at the given position or extrapolated value if position is out of range
109         * @see #setValue
110         */    
111        public float getValue(int position) {
112            return values[boundaryCondition.test(position, size())];
113        }
114    
115        /** 
116         * Set the signal value at the given position and deal with boundary conditions when position is out of range.
117         * Precondition : 0 &ge; <code>index</code> &le; size().
118         * @param   position     the position
119         * @param   signalValue  the value to set the signal to
120         * @see #getValue
121         */    
122        public void setValue(int position, float signalValue) {
123            values[boundaryCondition.test(position, size())] = signalValue;
124        }
125    
126        /** 
127         * Get a subset of the signal betwen the given start and end positions. The returned signal
128         * has the same size as the original signal but with values set to zero for out-of-range positions.
129         * @param   start   position where the subset starts
130         * @param   end     position where the subset ends
131         * @return          the subset of this signal in the given range
132         */
133        public SimpleSignal getSubset(int start, int end) {
134            SimpleSignal subset = new SimpleSignal(size());
135            for (int i = start; i < end; i++) subset.setValue(i, getValue(i));
136            return subset;
137        }
138    
139    
140        /** Get an estimate for the number of wavelet scales in this signal. 
141         * @return   the estimated number of wavelet scales contained in this signal, results are garanteed to be greater than zero
142         * @see Scale
143         */
144        public int estimateNumberOfScales() {
145            int maximumScale = 1 + (int) (log( 3.0d * ((double) size()) / 4.0d) / log(2.0d));
146            if (maximumScale >= 1) return maximumScale; else return 1;
147        }
148    
149        //-----------------factory methods--------------------------------------------------
150        
151        /** 
152         * Get a dicrete approximation of the Dirac delta function. There is an issue
153         * with signal having a size which is even : The middle is not defined !
154         * The signal has the value of one at the first position and zero everywhere else and
155         * periodic boundary conditions.   
156         * @return a dirac function signal 
157         */
158        public static SimpleSignal getDirac(int size) {
159            SimpleSignal dirac = new SimpleSignal(size);
160            dirac.setValue(0, 1f);
161            dirac.setBoundaryCondition(BoundaryCondition.Periodic);
162            return dirac;
163        }
164        
165        /** 
166         * Get a unit amplitude chirp signal with the given start and end frequencies.
167         * @param   size            the size of the signal
168         * @param   startFrequency  the starting frequency (in units of sample size)
169         * @param   endFrequency    the ending frequency (in units of sample size)
170         * @return                  the chirp
171         */
172        public static SimpleSignal getChirp(int size, float startFrequency, float endFrequency) {
173            SimpleSignal dirac = new SimpleSignal(size);
174            dirac.setValue(0, 1f);
175            return dirac;
176        }
177       
178        /** 
179         * Get a unit variance gaussian noise signal.
180         * @param   size            the size of the signal
181         * @return                  the gaussian noise signal
182         */
183        public static SimpleSignal getGaussianNoise(int size) {
184            SimpleSignal noise = new SimpleSignal(size);
185            for (int i = 0; i < size; i++) noise.setValue(i, getGaussianDeviate());
186            return noise;
187        }
188        
189        static int iset = 0;
190        static float gset;
191        
192        /** 
193         * Get a normally distributed gaussian deviate with zero mean and unit variance. 
194         * (see function gasdev(long *idnum) in Press (1992) Numerical Recipes in C, 2nd ed. page 289)
195         * @return a number roughly between -1 and +1
196         */
197        public static float getGaussianDeviate() {
198            if (iset == 0) {
199                float v1, v2, rsq;
200                do {  //TODO: use better random number generator like ran1(idnum) from Num.Rec.C
201                    v1 = 2f * (float) random() - 1f;  // Using java.lang.Math random number generator
202                    v2 = 2f * (float) random() - 1f;
203                    rsq = v1*v1 + v2*v2;
204                } while (rsq >= 1f || rsq == 0f);
205                float fac = (float) sqrt(-2d * log(rsq)/rsq);
206                gset = v1 * fac;
207                iset = 1;
208                return v2 * fac;
209            } else {
210                iset = 0;
211                return gset;
212            }
213        }
214    
215        /** 
216         * Get the B3-spline scaling funtion of the a-trous (with Mirror boundary)
217         * @return a scaling function
218         */
219        public static SimpleSignal getScalingFunction(int size, float scale) {
220            SimpleSignal scalingFunction = new SimpleSignal(size);
221            for (int i = 0; i < size; i++) {
222                float x = scale * i;
223                scalingFunction.setValue(i, (      aCube(x - 2f) 
224                                            - 4f * aCube(x - 1f) 
225                                            + 6f * aCube(x) 
226                                            - 4f * aCube(x + 1f)
227                                            +      aCube(x + 2f))/12f);
228            }
229                
230            return scalingFunction;
231        }
232    
233        static float aCube(float x) {
234            float absX = abs(x);
235            return absX * absX * absX;
236        }
237    
238        /** 
239         * Get the B3-spline wavelet funtion of the a-trous (with Mirror boundary)
240         * @return a wavelet function
241         */
242        public static SimpleSignal getWavelet(int size, float scale) {
243            return getScalingFunction(size, 2f * scale).multiply(2f).subtract(getScalingFunction(size, scale));
244        }
245    
246        //-----------------------simple math---------------------------------------------------
247    
248        /** 
249         * Add the given signal to this signals and return result in this instance.
250         * Precondition: signal must be non-null and equal in size to this signal
251         * @param signal the signal to be added
252         */
253        public SimpleSignal add(SimpleSignal signal) {
254            add(this, signal, this);
255            return this;
256        }
257    
258        /** 
259         * Subtract the given signal from this signal and return result in this instance.
260         * Precondition: signal must be non-null and equal in size to this signal
261         * @param signal the signal to be subtracted
262         */
263        public SimpleSignal subtract(SimpleSignal signal) {
264            subtract(this, signal, this);
265            return this;
266        }
267    
268        /** 
269         * Multiply by the given signal and return result in this instance.
270         * Precondition: signal must be non-null and equal in size to this signal
271         * @param signal the signal to be multiplied
272         */
273        public SimpleSignal multiply(SimpleSignal signal) {
274            multiply(this, signal, this);
275            return this;
276        }
277        
278        /** 
279         * Multiply by the given multipler and return result in this instance.
280         * Precondition: signal must be non-null and equal in size to this signal
281         * @param multiplier the multiplier
282         */
283        public SimpleSignal multiply(float multiplier) {
284            multiply(this, multiplier, this);
285            return this;
286        }
287        /** 
288         * Square the given signal and return result in this instance.
289         * Precondition: signal must be non-null and equal in size to this signal
290         */
291        public SimpleSignal square() {
292            square(this, this);
293            return this;
294        }
295    
296        /** 
297         * Return in 'c' the addition of 'a' and 'b'. 
298         * @param signalA input
299         * @param signalB input
300         * @param signalC output equals signalA + signalB
301         */
302        public static void add(SimpleSignal signalA, SimpleSignal signalB, SimpleSignal signalC) {
303            for (int i = 0; i < signalA.size(); i++) signalC.setValue(i, signalA.getValue(i) + signalB.getValue(i));
304        }
305    
306        /** 
307         * Return in 'c' the subtraction of 'a' minus 'b'. 
308         * @param signalA input
309         * @param signalB input
310         * @param signalC output equals signalA - signalB
311         */
312        public static void subtract(SimpleSignal signalA, SimpleSignal signalB, SimpleSignal signalC) {
313            for (int i = 0; i < signalA.size(); i++) signalC.setValue(i, signalA.getValue(i) - signalB.getValue(i));
314        }
315    
316        /** 
317         * Return in 'c' the multiplication of 'a' and 'b'. 
318         * @param signalA input
319         * @param signalB input
320         * @param signalC output equals signalA * signalB
321         */
322        public static void multiply(SimpleSignal signalA, SimpleSignal signalB, SimpleSignal signalC) {
323            for (int i = 0; i < signalA.size(); i++) signalC.setValue(i, signalA.getValue(i) * signalB.getValue(i));
324        }
325    
326        /** 
327         * Return in 'B' the multiplication of 'a' and multiplier. 
328         * @param signalA input
329         * @param multiplier input
330         * @param signalB output equals signalA * multiplier
331         */
332        public static void multiply(SimpleSignal signalA, float multiplier, SimpleSignal signalB) {
333            for (int i = 0; i < signalA.size(); i++) signalB.setValue(i, signalA.getValue(i) * multiplier);
334        }
335    
336        /** 
337         * Return in 'b' the square of 'a' and 'b'. 
338         * @param signalA input
339         * @param signalB output equals signalA*signalA
340         */
341        public static void square(SimpleSignal signalA, SimpleSignal signalB) {
342            for (int i = 0; i < signalA.size(); i++) signalB.setValue(i, signalA.getValue(i) * signalA.getValue(i));
343        }
344       
345    }