001/* ***** BEGIN LICENSE BLOCK *****
002 * Version: MPL 1.1/GPL 2.0/LGPL 2.1
003 *
004 * The contents of this file are subject to the Mozilla Public License Version
005 * 1.1 (the "License"); you may not use this file except in compliance with
006 * the License. You may obtain a copy of the License at
007 * http://www.mozilla.org/MPL/
008 *
009 * Software distributed under the License is distributed on an "AS IS" basis,
010 * WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License
011 * for the specific language governing rights and limitations under the
012 * License.
013 *
014 * The Original Code is JTransforms.
015 *
016 * The Initial Developer of the Original Code is
017 * Piotr Wendykier, Emory University.
018 * Portions created by the Initial Developer are Copyright (C) 2007-2009
019 * the Initial Developer. All Rights Reserved.
020 *
021 * Alternatively, the contents of this file may be used under the terms of
022 * either the GNU General Public License Version 2 or later (the "GPL"), or
023 * the GNU Lesser General Public License Version 2.1 or later (the "LGPL"),
024 * in which case the provisions of the GPL or the LGPL are applicable instead
025 * of those above. If you wish to allow use of your version of this file only
026 * under the terms of either the GPL or the LGPL, and not to allow others to
027 * use your version of this file under the terms of the MPL, indicate your
028 * decision by deleting the provisions above and replace them with the notice
029 * and other provisions required by the GPL or the LGPL. If you do not delete
030 * the provisions above, a recipient may use your version of this file under
031 * the terms of any one of the MPL, the GPL or the LGPL.
032 *
033 * ***** END LICENSE BLOCK ***** */
034
035package edu.emory.mathcs.jtransforms.dst;
036
037import java.util.concurrent.Future;
038
039import edu.emory.mathcs.jtransforms.dct.FloatDCT_1D;
040import edu.emory.mathcs.utils.ConcurrencyUtils;
041
042/**
043 * Computes 1D Discrete Sine Transform (DST) of single precision data. The size
044 * of data can be an arbitrary number. It uses DCT algorithm. This is a parallel
045 * implementation optimized for SMP systems.
046 * 
047 * @author Piotr Wendykier (piotr.wendykier@gmail.com)
048 * 
049 */
050public class FloatDST_1D {
051
052    private int n;
053
054    private FloatDCT_1D dct;
055
056    /**
057     * Creates new instance of FloatDST_1D.
058     * 
059     * @param n
060     *            size of data
061     */
062    public FloatDST_1D(int n) {
063        this.n = n;
064        dct = new FloatDCT_1D(n);
065    }
066
067    /**
068     * Computes 1D forward DST (DST-II) leaving the result in <code>a</code>.
069     * 
070     * @param a
071     *            data to transform
072     * @param scale
073     *            if true then scaling is performed
074     */
075    public void forward(float[] a, boolean scale) {
076        forward(a, 0, scale);
077    }
078
079    /**
080     * Computes 1D forward DST (DST-II) leaving the result in <code>a</code>.
081     * 
082     * @param a
083     *            data to transform
084     * @param offa
085     *            index of the first element in array <code>a</code>
086     * @param scale
087     *            if true then scaling is performed
088     */
089    public void forward(final float[] a, final int offa, boolean scale) {
090        if (n == 1)
091            return;
092        float tmp;
093        int nd2 = n / 2;
094        int startIdx = 1 + offa;
095        int stopIdx = offa + n;
096        for (int i = startIdx; i < stopIdx; i += 2) {
097            a[i] = -a[i];
098        }
099        dct.forward(a, offa, scale);
100        int nthreads = ConcurrencyUtils.getNumberOfThreads();
101        if ((nthreads > 1) && (nd2 > ConcurrencyUtils.getThreadsBeginN_1D_FFT_2Threads())) {
102            nthreads = 2;
103            final int k = nd2 / nthreads;
104            Future<?>[] futures = new Future[nthreads];
105            for (int j = 0; j < nthreads; j++) {
106                final int firstIdx = j * k;
107                final int lastIdx = (j == (nthreads - 1)) ? nd2 : firstIdx + k;
108                futures[j] = ConcurrencyUtils.submit(new Runnable() {
109                    public void run() {
110                        float tmp;
111                        int idx0 = offa + n - 1;
112                        int idx1;
113                        int idx2;
114                        for (int i = firstIdx; i < lastIdx; i++) {
115                            idx2 = offa + i;
116                            tmp = a[idx2];
117                            idx1 = idx0 - i;
118                            a[idx2] = a[idx1];
119                            a[idx1] = tmp;
120                        }
121                    }
122                });
123            }
124            ConcurrencyUtils.waitForCompletion(futures);
125        } else {
126            int idx0 = offa + n - 1;
127            int idx1;
128            int idx2;
129            for (int i = 0; i < nd2; i++) {
130                idx2 = offa + i;
131                tmp = a[idx2];
132                idx1 = idx0 - i;
133                a[idx2] = a[idx1];
134                a[idx1] = tmp;
135            }
136        }
137    }
138
139    /**
140     * Computes 1D inverse DST (DST-III) leaving the result in <code>a</code>.
141     * 
142     * @param a
143     *            data to transform
144     * @param scale
145     *            if true then scaling is performed
146     */
147    public void inverse(float[] a, boolean scale) {
148        inverse(a, 0, scale);
149    }
150
151    /**
152     * Computes 1D inverse DST (DST-III) leaving the result in <code>a</code>.
153     * 
154     * @param a
155     *            data to transform
156     * @param offa
157     *            index of the first element in array <code>a</code>
158     * @param scale
159     *            if true then scaling is performed
160     */
161    public void inverse(final float[] a, final int offa, boolean scale) {
162        if (n == 1)
163            return;
164        float tmp;
165        int nd2 = n / 2;
166        int nthreads = ConcurrencyUtils.getNumberOfThreads();
167        if ((nthreads > 1) && (nd2 > ConcurrencyUtils.getThreadsBeginN_1D_FFT_2Threads())) {
168            nthreads = 2;
169            final int k = nd2 / nthreads;
170            Future<?>[] futures = new Future[nthreads];
171            for (int j = 0; j < nthreads; j++) {
172                final int firstIdx = j * k;
173                final int lastIdx = (j == (nthreads - 1)) ? nd2 : firstIdx + k;
174                futures[j] = ConcurrencyUtils.submit(new Runnable() {
175                    public void run() {
176                        float tmp;
177                        int idx0 = offa + n - 1;
178                        int idx1, idx2;
179                        for (int i = firstIdx; i < lastIdx; i++) {
180                            idx2 = offa + i;
181                            tmp = a[idx2];
182                            idx1 = idx0 - i;
183                            a[idx2] = a[idx1];
184                            a[idx1] = tmp;
185                        }
186                    }
187                });
188            }
189            ConcurrencyUtils.waitForCompletion(futures);
190        } else {
191            int idx0 = offa + n - 1;
192            for (int i = 0; i < nd2; i++) {
193                tmp = a[offa + i];
194                a[offa + i] = a[idx0 - i];
195                a[idx0 - i] = tmp;
196            }
197        }
198        dct.inverse(a, offa, scale);
199        int startidx = 1 + offa;
200        int stopidx = offa + n;
201        for (int i = startidx; i < stopidx; i += 2) {
202            a[i] = -a[i];
203        }
204    }
205}