1/* 2 * Licensed to the Apache Software Foundation (ASF) under one or more 3 * contributor license agreements. See the NOTICE file distributed with 4 * this work for additional information regarding copyright ownership. 5 * The ASF licenses this file to You under the Apache License, Version 2.0 6 * (the "License"); you may not use this file except in compliance with 7 * the License. You may obtain a copy of the License at 8 * 9 * http://www.apache.org/licenses/LICENSE-2.0 10 * 11 * Unless required by applicable law or agreed to in writing, software 12 * distributed under the License is distributed on an "AS IS" BASIS, 13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14 * See the License for the specific language governing permissions and 15 * limitations under the License. 16 */ 17package org.apache.commons.math.transform; 18 19import org.apache.commons.math.FunctionEvaluationException; 20import org.apache.commons.math.MathRuntimeException; 21import org.apache.commons.math.analysis.UnivariateRealFunction; 22import org.apache.commons.math.complex.Complex; 23import org.apache.commons.math.exception.util.LocalizedFormats; 24import org.apache.commons.math.util.FastMath; 25 26/** 27 * Implements the <a href="http://documents.wolfram.com/v5/Add-onsLinks/ 28 * StandardPackages/LinearAlgebra/FourierTrig.html">Fast Sine Transform</a> 29 * for transformation of one-dimensional data sets. For reference, see 30 * <b>Fast Fourier Transforms</b>, ISBN 0849371635, chapter 3. 31 * <p> 32 * FST is its own inverse, up to a multiplier depending on conventions. 33 * The equations are listed in the comments of the corresponding methods.</p> 34 * <p> 35 * Similar to FFT, we also require the length of data set to be power of 2. 36 * In addition, the first element must be 0 and it's enforced in function 37 * transformation after sampling.</p> 38 * <p>As of version 2.0 this no longer implements Serializable</p> 39 * 40 * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 févr. 2011) $ 41 * @since 1.2 42 */ 43public class FastSineTransformer implements RealTransformer { 44 45 /** 46 * Construct a default transformer. 47 */ 48 public FastSineTransformer() { 49 super(); 50 } 51 52 /** 53 * Transform the given real data set. 54 * <p> 55 * The formula is F<sub>n</sub> = ∑<sub>k=0</sub><sup>N-1</sup> f<sub>k</sub> sin(π nk/N) 56 * </p> 57 * 58 * @param f the real data array to be transformed 59 * @return the real transformed array 60 * @throws IllegalArgumentException if any parameters are invalid 61 */ 62 public double[] transform(double f[]) 63 throws IllegalArgumentException { 64 return fst(f); 65 } 66 67 /** 68 * Transform the given real function, sampled on the given interval. 69 * <p> 70 * The formula is F<sub>n</sub> = ∑<sub>k=0</sub><sup>N-1</sup> f<sub>k</sub> sin(π nk/N) 71 * </p> 72 * 73 * @param f the function to be sampled and transformed 74 * @param min the lower bound for the interval 75 * @param max the upper bound for the interval 76 * @param n the number of sample points 77 * @return the real transformed array 78 * @throws FunctionEvaluationException if function cannot be evaluated 79 * at some point 80 * @throws IllegalArgumentException if any parameters are invalid 81 */ 82 public double[] transform(UnivariateRealFunction f, 83 double min, double max, int n) 84 throws FunctionEvaluationException, IllegalArgumentException { 85 86 double data[] = FastFourierTransformer.sample(f, min, max, n); 87 data[0] = 0.0; 88 return fst(data); 89 } 90 91 /** 92 * Transform the given real data set. 93 * <p> 94 * The formula is F<sub>n</sub> = √(2/N) ∑<sub>k=0</sub><sup>N-1</sup> f<sub>k</sub> sin(π nk/N) 95 * </p> 96 * 97 * @param f the real data array to be transformed 98 * @return the real transformed array 99 * @throws IllegalArgumentException if any parameters are invalid 100 */ 101 public double[] transform2(double f[]) throws IllegalArgumentException { 102 103 double scaling_coefficient = FastMath.sqrt(2.0 / f.length); 104 return FastFourierTransformer.scaleArray(fst(f), scaling_coefficient); 105 } 106 107 /** 108 * Transform the given real function, sampled on the given interval. 109 * <p> 110 * The formula is F<sub>n</sub> = √(2/N) ∑<sub>k=0</sub><sup>N-1</sup> f<sub>k</sub> sin(π nk/N) 111 * </p> 112 * 113 * @param f the function to be sampled and transformed 114 * @param min the lower bound for the interval 115 * @param max the upper bound for the interval 116 * @param n the number of sample points 117 * @return the real transformed array 118 * @throws FunctionEvaluationException if function cannot be evaluated 119 * at some point 120 * @throws IllegalArgumentException if any parameters are invalid 121 */ 122 public double[] transform2( 123 UnivariateRealFunction f, double min, double max, int n) 124 throws FunctionEvaluationException, IllegalArgumentException { 125 126 double data[] = FastFourierTransformer.sample(f, min, max, n); 127 data[0] = 0.0; 128 double scaling_coefficient = FastMath.sqrt(2.0 / n); 129 return FastFourierTransformer.scaleArray(fst(data), scaling_coefficient); 130 } 131 132 /** 133 * Inversely transform the given real data set. 134 * <p> 135 * The formula is f<sub>k</sub> = (2/N) ∑<sub>n=0</sub><sup>N-1</sup> F<sub>n</sub> sin(π nk/N) 136 * </p> 137 * 138 * @param f the real data array to be inversely transformed 139 * @return the real inversely transformed array 140 * @throws IllegalArgumentException if any parameters are invalid 141 */ 142 public double[] inversetransform(double f[]) throws IllegalArgumentException { 143 144 double scaling_coefficient = 2.0 / f.length; 145 return FastFourierTransformer.scaleArray(fst(f), scaling_coefficient); 146 } 147 148 /** 149 * Inversely transform the given real function, sampled on the given interval. 150 * <p> 151 * The formula is f<sub>k</sub> = (2/N) ∑<sub>n=0</sub><sup>N-1</sup> F<sub>n</sub> sin(π nk/N) 152 * </p> 153 * 154 * @param f the function to be sampled and inversely transformed 155 * @param min the lower bound for the interval 156 * @param max the upper bound for the interval 157 * @param n the number of sample points 158 * @return the real inversely transformed array 159 * @throws FunctionEvaluationException if function cannot be evaluated at some point 160 * @throws IllegalArgumentException if any parameters are invalid 161 */ 162 public double[] inversetransform(UnivariateRealFunction f, double min, double max, int n) 163 throws FunctionEvaluationException, IllegalArgumentException { 164 165 double data[] = FastFourierTransformer.sample(f, min, max, n); 166 data[0] = 0.0; 167 double scaling_coefficient = 2.0 / n; 168 return FastFourierTransformer.scaleArray(fst(data), scaling_coefficient); 169 } 170 171 /** 172 * Inversely transform the given real data set. 173 * <p> 174 * The formula is f<sub>k</sub> = √(2/N) ∑<sub>n=0</sub><sup>N-1</sup> F<sub>n</sub> sin(π nk/N) 175 * </p> 176 * 177 * @param f the real data array to be inversely transformed 178 * @return the real inversely transformed array 179 * @throws IllegalArgumentException if any parameters are invalid 180 */ 181 public double[] inversetransform2(double f[]) throws IllegalArgumentException { 182 183 return transform2(f); 184 } 185 186 /** 187 * Inversely transform the given real function, sampled on the given interval. 188 * <p> 189 * The formula is f<sub>k</sub> = √(2/N) ∑<sub>n=0</sub><sup>N-1</sup> F<sub>n</sub> sin(π nk/N) 190 * </p> 191 * 192 * @param f the function to be sampled and inversely transformed 193 * @param min the lower bound for the interval 194 * @param max the upper bound for the interval 195 * @param n the number of sample points 196 * @return the real inversely transformed array 197 * @throws FunctionEvaluationException if function cannot be evaluated at some point 198 * @throws IllegalArgumentException if any parameters are invalid 199 */ 200 public double[] inversetransform2(UnivariateRealFunction f, double min, double max, int n) 201 throws FunctionEvaluationException, IllegalArgumentException { 202 203 return transform2(f, min, max, n); 204 } 205 206 /** 207 * Perform the FST algorithm (including inverse). 208 * 209 * @param f the real data array to be transformed 210 * @return the real transformed array 211 * @throws IllegalArgumentException if any parameters are invalid 212 */ 213 protected double[] fst(double f[]) throws IllegalArgumentException { 214 215 final double transformed[] = new double[f.length]; 216 217 FastFourierTransformer.verifyDataSet(f); 218 if (f[0] != 0.0) { 219 throw MathRuntimeException.createIllegalArgumentException( 220 LocalizedFormats.FIRST_ELEMENT_NOT_ZERO, 221 f[0]); 222 } 223 final int n = f.length; 224 if (n == 1) { // trivial case 225 transformed[0] = 0.0; 226 return transformed; 227 } 228 229 // construct a new array and perform FFT on it 230 final double[] x = new double[n]; 231 x[0] = 0.0; 232 x[n >> 1] = 2.0 * f[n >> 1]; 233 for (int i = 1; i < (n >> 1); i++) { 234 final double a = FastMath.sin(i * FastMath.PI / n) * (f[i] + f[n-i]); 235 final double b = 0.5 * (f[i] - f[n-i]); 236 x[i] = a + b; 237 x[n - i] = a - b; 238 } 239 FastFourierTransformer transformer = new FastFourierTransformer(); 240 Complex y[] = transformer.transform(x); 241 242 // reconstruct the FST result for the original array 243 transformed[0] = 0.0; 244 transformed[1] = 0.5 * y[0].getReal(); 245 for (int i = 1; i < (n >> 1); i++) { 246 transformed[2 * i] = -y[i].getImaginary(); 247 transformed[2 * i + 1] = y[i].getReal() + transformed[2 * i - 1]; 248 } 249 250 return transformed; 251 } 252} 253