001 /* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017 package org.apache.commons.math.transform; 018 019 import org.apache.commons.math.FunctionEvaluationException; 020 import org.apache.commons.math.MathRuntimeException; 021 import org.apache.commons.math.analysis.UnivariateRealFunction; 022 import org.apache.commons.math.complex.Complex; 023 import org.apache.commons.math.exception.util.LocalizedFormats; 024 import org.apache.commons.math.util.FastMath; 025 026 /** 027 * Implements the <a href="http://documents.wolfram.com/v5/Add-onsLinks/ 028 * StandardPackages/LinearAlgebra/FourierTrig.html">Fast Sine Transform</a> 029 * for transformation of one-dimensional data sets. For reference, see 030 * <b>Fast Fourier Transforms</b>, ISBN 0849371635, chapter 3. 031 * <p> 032 * FST is its own inverse, up to a multiplier depending on conventions. 033 * The equations are listed in the comments of the corresponding methods.</p> 034 * <p> 035 * Similar to FFT, we also require the length of data set to be power of 2. 036 * In addition, the first element must be 0 and it's enforced in function 037 * transformation after sampling.</p> 038 * <p>As of version 2.0 this no longer implements Serializable</p> 039 * 040 * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 f??vr. 2011) $ 041 * @since 1.2 042 */ 043 public class FastSineTransformer implements RealTransformer { 044 045 /** 046 * Construct a default transformer. 047 */ 048 public FastSineTransformer() { 049 super(); 050 } 051 052 /** 053 * Transform the given real data set. 054 * <p> 055 * The formula is F<sub>n</sub> = ∑<sub>k=0</sub><sup>N-1</sup> f<sub>k</sub> sin(π nk/N) 056 * </p> 057 * 058 * @param f the real data array to be transformed 059 * @return the real transformed array 060 * @throws IllegalArgumentException if any parameters are invalid 061 */ 062 public double[] transform(double f[]) 063 throws IllegalArgumentException { 064 return fst(f); 065 } 066 067 /** 068 * Transform the given real function, sampled on the given interval. 069 * <p> 070 * The formula is F<sub>n</sub> = ∑<sub>k=0</sub><sup>N-1</sup> f<sub>k</sub> sin(π nk/N) 071 * </p> 072 * 073 * @param f the function to be sampled and transformed 074 * @param min the lower bound for the interval 075 * @param max the upper bound for the interval 076 * @param n the number of sample points 077 * @return the real transformed array 078 * @throws FunctionEvaluationException if function cannot be evaluated 079 * at some point 080 * @throws IllegalArgumentException if any parameters are invalid 081 */ 082 public double[] transform(UnivariateRealFunction f, 083 double min, double max, int n) 084 throws FunctionEvaluationException, IllegalArgumentException { 085 086 double data[] = FastFourierTransformer.sample(f, min, max, n); 087 data[0] = 0.0; 088 return fst(data); 089 } 090 091 /** 092 * Transform the given real data set. 093 * <p> 094 * The formula is F<sub>n</sub> = √(2/N) ∑<sub>k=0</sub><sup>N-1</sup> f<sub>k</sub> sin(π nk/N) 095 * </p> 096 * 097 * @param f the real data array to be transformed 098 * @return the real transformed array 099 * @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 }