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 */ 017package org.apache.commons.math4.transform; 018 019import java.util.function.UnaryOperator; 020import java.util.function.DoubleUnaryOperator; 021 022import org.apache.commons.numbers.complex.Complex; 023import org.apache.commons.numbers.core.ArithmeticUtils; 024 025/** 026 * Implements the Fast Sine Transform for transformation of one-dimensional real 027 * data sets. For reference, see James S. Walker, <em>Fast Fourier 028 * Transforms</em>, chapter 3 (ISBN 0849371635). 029 * <p> 030 * There are several variants of the discrete sine transform. The present 031 * implementation corresponds to DST-I, with various normalization conventions, 032 * which are specified by the parameter {@link Norm}. 033 * <strong>It should be noted that regardless to the convention, the first 034 * element of the dataset to be transformed must be zero.</strong> 035 * <p> 036 * DST-I is equivalent to DFT of an <em>odd extension</em> of the data series. 037 * More precisely, if x<sub>0</sub>, …, x<sub>N-1</sub> is the data set 038 * to be sine transformed, the extended data set x<sub>0</sub><sup>#</sup>, 039 * …, x<sub>2N-1</sub><sup>#</sup> is defined as follows 040 * <ul> 041 * <li>x<sub>0</sub><sup>#</sup> = x<sub>0</sub> = 0,</li> 042 * <li>x<sub>k</sub><sup>#</sup> = x<sub>k</sub> if 1 ≤ k < N,</li> 043 * <li>x<sub>N</sub><sup>#</sup> = 0,</li> 044 * <li>x<sub>k</sub><sup>#</sup> = -x<sub>2N-k</sub> if N + 1 ≤ k < 045 * 2N.</li> 046 * </ul> 047 * <p> 048 * Then, the standard DST-I y<sub>0</sub>, …, y<sub>N-1</sub> of the real 049 * data set x<sub>0</sub>, …, x<sub>N-1</sub> is equal to <em>half</em> 050 * of i (the pure imaginary number) times the N first elements of the DFT of the 051 * extended data set x<sub>0</sub><sup>#</sup>, …, 052 * x<sub>2N-1</sub><sup>#</sup> <br> 053 * y<sub>n</sub> = (i / 2) ∑<sub>k=0</sub><sup>2N-1</sup> 054 * x<sub>k</sub><sup>#</sup> exp[-2πi nk / (2N)] 055 * k = 0, …, N-1. 056 * <p> 057 * The present implementation of the discrete sine transform as a fast sine 058 * transform requires the length of the data to be a power of two. Besides, 059 * it implicitly assumes that the sampled function is odd. In particular, the 060 * first element of the data set must be 0, which is enforced in 061 * {@link #apply(DoubleUnaryOperator, double, double, int)}, 062 * after sampling. 063 */ 064public class FastSineTransform implements RealTransform { 065 /** Operation to be performed. */ 066 private final UnaryOperator<double[]> op; 067 068 /** 069 * @param normalization Normalization to be applied to the transformed data. 070 * @param inverse Whether to perform the inverse transform. 071 */ 072 public FastSineTransform(final Norm normalization, 073 final boolean inverse) { 074 op = create(normalization, inverse); 075 } 076 077 /** 078 * @param normalization Normalization to be applied to the 079 * transformed data. 080 */ 081 public FastSineTransform(final Norm normalization) { 082 this(normalization, false); 083 } 084 085 /** 086 * {@inheritDoc} 087 * 088 * The first element of the specified data set is required to be {@code 0}. 089 * 090 * @throws IllegalArgumentException if the length of the data array is 091 * not a power of two, or the first element of the data array is not zero. 092 */ 093 @Override 094 public double[] apply(final double[] f) { 095 return op.apply(f); 096 } 097 098 /** 099 * {@inheritDoc} 100 * 101 * The implementation enforces {@code f(x) = 0} at {@code x = 0}. 102 * 103 * @throws IllegalArgumentException if the number of sample points is not a 104 * power of two, if the lower bound is greater than, or equal to the upper bound, 105 * if the number of sample points is negative. 106 */ 107 @Override 108 public double[] apply(final DoubleUnaryOperator f, 109 final double min, 110 final double max, 111 final int n) { 112 final double[] data = TransformUtils.sample(f, min, max, n); 113 data[0] = 0; 114 return apply(data); 115 } 116 117 /** 118 * Perform the FST algorithm (including inverse). 119 * The first element of the data set is required to be {@code 0}. 120 * 121 * @param f Data array to be transformed. 122 * @return the transformed array. 123 * @throws IllegalArgumentException if the length of the data array is 124 * not a power of two, or the first element of the data array is not zero. 125 */ 126 private double[] fst(double[] f) { 127 if (!ArithmeticUtils.isPowerOfTwo(f.length)) { 128 throw new TransformException(TransformException.NOT_POWER_OF_TWO, 129 f.length); 130 } 131 if (f[0] != 0) { 132 throw new TransformException(TransformException.FIRST_ELEMENT_NOT_ZERO, 133 f[0]); 134 } 135 136 final double[] transformed = new double[f.length]; 137 final int n = f.length; 138 if (n == 1) { 139 transformed[0] = 0; 140 return transformed; 141 } 142 143 // construct a new array and perform FFT on it 144 final double[] x = new double[n]; 145 x[0] = 0; 146 final int nShifted = n >> 1; 147 x[nShifted] = 2 * f[nShifted]; 148 final double piOverN = Math.PI / n; 149 for (int i = 1; i < nShifted; i++) { 150 final int nMi = n - i; 151 final double fi = f[i]; 152 final double fnMi = f[nMi]; 153 final double a = Math.sin(i * piOverN) * (fi + fnMi); 154 final double b = 0.5 * (fi - fnMi); 155 x[i] = a + b; 156 x[nMi] = a - b; 157 } 158 159 final FastFourierTransform transform = new FastFourierTransform(FastFourierTransform.Norm.STD); 160 final Complex[] y = transform.apply(x); 161 162 // reconstruct the FST result for the original array 163 transformed[0] = 0; 164 transformed[1] = 0.5 * y[0].getReal(); 165 for (int i = 1; i < nShifted; i++) { 166 final int i2 = 2 * i; 167 transformed[i2] = -y[i].getImaginary(); 168 transformed[i2 + 1] = y[i].getReal() + transformed[i2 - 1]; 169 } 170 171 return transformed; 172 } 173 174 /** 175 * Factory method. 176 * 177 * @param normalization Normalization to be applied to the 178 * transformed data. 179 * @param inverse Whether to perform the inverse transform. 180 * @return the transform operator. 181 */ 182 private UnaryOperator<double[]> create(final Norm normalization, 183 final boolean inverse) { 184 if (inverse) { 185 return normalization == Norm.ORTHO ? 186 f -> TransformUtils.scaleInPlace(fst(f), Math.sqrt(2d / f.length)) : 187 f -> TransformUtils.scaleInPlace(fst(f), 2d / f.length); 188 } else { 189 return normalization == Norm.ORTHO ? 190 f -> TransformUtils.scaleInPlace(fst(f), Math.sqrt(2d / f.length)) : 191 f -> fst(f); 192 } 193 } 194 195 /** 196 * Normalization types. 197 */ 198 public enum Norm { 199 /** 200 * Should be passed to the constructor of {@link FastSineTransform} to 201 * use the <em>standard</em> normalization convention. The standard DST-I 202 * normalization convention is defined as follows 203 * <ul> 204 * <li>forward transform: y<sub>n</sub> = ∑<sub>k=0</sub><sup>N-1</sup> 205 * x<sub>k</sub> sin(π nk / N),</li> 206 * <li>inverse transform: x<sub>k</sub> = (2 / N) 207 * ∑<sub>n=0</sub><sup>N-1</sup> y<sub>n</sub> sin(π nk / N),</li> 208 * </ul> 209 * where N is the size of the data sample, and x<sub>0</sub> = 0. 210 */ 211 STD, 212 213 /** 214 * Should be passed to the constructor of {@link FastSineTransform} to 215 * use the <em>orthogonal</em> normalization convention. The orthogonal 216 * DCT-I normalization convention is defined as follows 217 * <ul> 218 * <li>Forward transform: y<sub>n</sub> = √(2 / N) 219 * ∑<sub>k=0</sub><sup>N-1</sup> x<sub>k</sub> sin(π nk / N),</li> 220 * <li>Inverse transform: x<sub>k</sub> = √(2 / N) 221 * ∑<sub>n=0</sub><sup>N-1</sup> y<sub>n</sub> sin(π nk / N),</li> 222 * </ul> 223 * which makes the transform orthogonal. N is the size of the data sample, 224 * and x<sub>0</sub> = 0. 225 */ 226 ORTHO 227 } 228}