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.analysis.*; 020 import org.apache.commons.math.complex.*; 021 import org.apache.commons.math.FunctionEvaluationException; 022 import org.apache.commons.math.MathRuntimeException; 023 024 /** 025 * Implements the <a href="http://documents.wolfram.com/v5/Add-onsLinks/ 026 * StandardPackages/LinearAlgebra/FourierTrig.html">Fast Sine Transform</a> 027 * for transformation of one-dimensional data sets. For reference, see 028 * <b>Fast Fourier Transforms</b>, ISBN 0849371635, chapter 3. 029 * <p> 030 * FST is its own inverse, up to a multiplier depending on conventions. 031 * The equations are listed in the comments of the corresponding methods.</p> 032 * <p> 033 * Similar to FFT, we also require the length of data set to be power of 2. 034 * In addition, the first element must be 0 and it's enforced in function 035 * transformation after sampling.</p> 036 * <p>As of version 2.0 this no longer implements Serializable</p> 037 * 038 * @version $Revision: 780541 $ $Date: 2009-05-31 20:47:02 -0400 (Sun, 31 May 2009) $ 039 * @since 1.2 040 */ 041 public class FastSineTransformer implements RealTransformer { 042 043 /** 044 * Construct a default transformer. 045 */ 046 public FastSineTransformer() { 047 super(); 048 } 049 050 /** 051 * Transform the given real data set. 052 * <p> 053 * The formula is F<sub>n</sub> = ∑<sub>k=0</sub><sup>N-1</sup> f<sub>k</sub> sin(π nk/N) 054 * </p> 055 * 056 * @param f the real data array to be transformed 057 * @return the real transformed array 058 * @throws IllegalArgumentException if any parameters are invalid 059 */ 060 public double[] transform(double f[]) 061 throws IllegalArgumentException { 062 return fst(f); 063 } 064 065 /** 066 * Transform the given real function, sampled on the given interval. 067 * <p> 068 * The formula is F<sub>n</sub> = ∑<sub>k=0</sub><sup>N-1</sup> f<sub>k</sub> sin(π nk/N) 069 * </p> 070 * 071 * @param f the function to be sampled and transformed 072 * @param min the lower bound for the interval 073 * @param max the upper bound for the interval 074 * @param n the number of sample points 075 * @return the real transformed array 076 * @throws FunctionEvaluationException if function cannot be evaluated 077 * at some point 078 * @throws IllegalArgumentException if any parameters are invalid 079 */ 080 public double[] transform(UnivariateRealFunction f, 081 double min, double max, int n) 082 throws FunctionEvaluationException, IllegalArgumentException { 083 084 double data[] = FastFourierTransformer.sample(f, min, max, n); 085 data[0] = 0.0; 086 return fst(data); 087 } 088 089 /** 090 * Transform the given real data set. 091 * <p> 092 * The formula is F<sub>n</sub> = √(2/N) ∑<sub>k=0</sub><sup>N-1</sup> f<sub>k</sub> sin(π nk/N) 093 * </p> 094 * 095 * @param f the real data array to be transformed 096 * @return the real transformed array 097 * @throws IllegalArgumentException if any parameters are invalid 098 */ 099 public double[] transform2(double f[]) throws IllegalArgumentException { 100 101 double scaling_coefficient = Math.sqrt(2.0 / f.length); 102 return FastFourierTransformer.scaleArray(fst(f), scaling_coefficient); 103 } 104 105 /** 106 * Transform the given real function, sampled on the given interval. 107 * <p> 108 * The formula is F<sub>n</sub> = √(2/N) ∑<sub>k=0</sub><sup>N-1</sup> f<sub>k</sub> sin(π nk/N) 109 * </p> 110 * 111 * @param f the function to be sampled and transformed 112 * @param min the lower bound for the interval 113 * @param max the upper bound for the interval 114 * @param n the number of sample points 115 * @return the real transformed array 116 * @throws FunctionEvaluationException if function cannot be evaluated 117 * at some point 118 * @throws IllegalArgumentException if any parameters are invalid 119 */ 120 public double[] transform2( 121 UnivariateRealFunction f, double min, double max, int n) 122 throws FunctionEvaluationException, IllegalArgumentException { 123 124 double data[] = FastFourierTransformer.sample(f, min, max, n); 125 data[0] = 0.0; 126 double scaling_coefficient = Math.sqrt(2.0 / n); 127 return FastFourierTransformer.scaleArray(fst(data), scaling_coefficient); 128 } 129 130 /** 131 * Inversely transform the given real data set. 132 * <p> 133 * The formula is f<sub>k</sub> = (2/N) ∑<sub>n=0</sub><sup>N-1</sup> F<sub>n</sub> sin(π nk/N) 134 * </p> 135 * 136 * @param f the real data array to be inversely transformed 137 * @return the real inversely transformed array 138 * @throws IllegalArgumentException if any parameters are invalid 139 */ 140 public double[] inversetransform(double f[]) throws IllegalArgumentException { 141 142 double scaling_coefficient = 2.0 / f.length; 143 return FastFourierTransformer.scaleArray(fst(f), scaling_coefficient); 144 } 145 146 /** 147 * Inversely transform the given real function, sampled on the given interval. 148 * <p> 149 * The formula is f<sub>k</sub> = (2/N) ∑<sub>n=0</sub><sup>N-1</sup> F<sub>n</sub> sin(π nk/N) 150 * </p> 151 * 152 * @param f the function to be sampled and inversely transformed 153 * @param min the lower bound for the interval 154 * @param max the upper bound for the interval 155 * @param n the number of sample points 156 * @return the real inversely transformed array 157 * @throws FunctionEvaluationException if function cannot be evaluated 158 * at some point 159 * @throws IllegalArgumentException if any parameters are invalid 160 */ 161 public double[] inversetransform(UnivariateRealFunction f, double min, double max, int n) 162 throws FunctionEvaluationException, IllegalArgumentException { 163 164 double data[] = FastFourierTransformer.sample(f, min, max, n); 165 data[0] = 0.0; 166 double scaling_coefficient = 2.0 / n; 167 return FastFourierTransformer.scaleArray(fst(data), scaling_coefficient); 168 } 169 170 /** 171 * Inversely transform the given real data set. 172 * <p> 173 * The formula is f<sub>k</sub> = √(2/N) ∑<sub>n=0</sub><sup>N-1</sup> F<sub>n</sub> sin(π nk/N) 174 * </p> 175 * 176 * @param f the real data array to be inversely transformed 177 * @return the real inversely transformed array 178 * @throws IllegalArgumentException if any parameters are invalid 179 */ 180 public double[] inversetransform2(double f[]) throws IllegalArgumentException { 181 182 return transform2(f); 183 } 184 185 /** 186 * Inversely transform the given real function, sampled on the given interval. 187 * <p> 188 * The formula is f<sub>k</sub> = √(2/N) ∑<sub>n=0</sub><sup>N-1</sup> F<sub>n</sub> sin(π nk/N) 189 * </p> 190 * 191 * @param f the function to be sampled and inversely transformed 192 * @param min the lower bound for the interval 193 * @param max the upper bound for the interval 194 * @param n the number of sample points 195 * @return the real inversely transformed array 196 * @throws FunctionEvaluationException if function cannot be evaluated 197 * 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 double A, B, x[], F[] = new double[f.length]; 216 217 FastFourierTransformer.verifyDataSet(f); 218 if (f[0] != 0.0) { 219 throw MathRuntimeException.createIllegalArgumentException( 220 "first element is not 0: {0}", 221 f[0]); 222 } 223 int N = f.length; 224 if (N == 1) { // trivial case 225 F[0] = 0.0; 226 return F; 227 } 228 229 // construct a new array and perform FFT on it 230 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 A = Math.sin(i * Math.PI / N) * (f[i] + f[N-i]); 235 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 F[0] = 0.0; 244 F[1] = 0.5 * y[0].getReal(); 245 for (int i = 1; i < (N >> 1); i++) { 246 F[2*i] = -y[i].getImaginary(); 247 F[2*i+1] = y[i].getReal() + F[2*i-1]; 248 } 249 250 return F; 251 } 252 }