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.DoubleDCT_1D;
040import edu.emory.mathcs.utils.ConcurrencyUtils;
041
042/**
043 * Computes 1D Discrete Sine Transform (DST) of double 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 DoubleDST_1D {
051
052    private int n;
053
054    private DoubleDCT_1D dct;
055
056    /**
057     * Creates new instance of DoubleDST_1D.
058     * 
059     * @param n
060     *            size of data
061     */
062    public DoubleDST_1D(int n) {
063        this.n = n;
064        dct = new DoubleDCT_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(double[] 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 double[] a, final int offa, boolean scale) {
090        if (n == 1)
091            return;
092        double 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                    @Override
110                                        public void run() {
111                        double tmp;
112                        int idx0 = offa + n - 1;
113                        int idx1;
114                        int idx2;
115                        for (int i = firstIdx; i < lastIdx; i++) {
116                            idx2 = offa + i;
117                            tmp = a[idx2];
118                            idx1 = idx0 - i;
119                            a[idx2] = a[idx1];
120                            a[idx1] = tmp;
121                        }
122                    }
123                });
124            }
125            ConcurrencyUtils.waitForCompletion(futures);
126        } else {
127            int idx0 = offa + n - 1;
128            int idx1;
129            int idx2;
130            for (int i = 0; i < nd2; i++) {
131                idx2 = offa + i;
132                tmp = a[idx2];
133                idx1 = idx0 - i;
134                a[idx2] = a[idx1];
135                a[idx1] = tmp;
136            }
137        }
138    }
139
140    /**
141     * Computes 1D inverse DST (DST-III) leaving the result in <code>a</code>.
142     * 
143     * @param a
144     *            data to transform
145     * @param scale
146     *            if true then scaling is performed
147     */
148    public void inverse(double[] a, boolean scale) {
149        inverse(a, 0, scale);
150    }
151
152    /**
153     * Computes 1D inverse DST (DST-III) leaving the result in <code>a</code>.
154     * 
155     * @param a
156     *            data to transform
157     * @param offa
158     *            index of the first element in array <code>a</code>
159     * @param scale
160     *            if true then scaling is performed
161     */
162    public void inverse(final double[] a, final int offa, boolean scale) {
163        if (n == 1)
164            return;
165        double tmp;
166        int nd2 = n / 2;
167        int nthreads = ConcurrencyUtils.getNumberOfThreads();
168        if ((nthreads > 1) && (nd2 > ConcurrencyUtils.getThreadsBeginN_1D_FFT_2Threads())) {
169            nthreads = 2;
170            final int k = nd2 / nthreads;
171            Future<?>[] futures = new Future[nthreads];
172            for (int j = 0; j < nthreads; j++) {
173                final int firstIdx = j * k;
174                final int lastIdx = (j == (nthreads - 1)) ? nd2 : firstIdx + k;
175                futures[j] = ConcurrencyUtils.submit(new Runnable() {
176                    @Override
177                                        public void run() {
178                        double tmp;
179                        int idx0 = offa + n - 1;
180                        int idx1, idx2;
181                        for (int i = firstIdx; i < lastIdx; i++) {
182                            idx2 = offa + i;
183                            tmp = a[idx2];
184                            idx1 = idx0 - i;
185                            a[idx2] = a[idx1];
186                            a[idx1] = tmp;
187                        }
188                    }
189                });
190            }
191            ConcurrencyUtils.waitForCompletion(futures);
192        } else {
193            int idx0 = offa + n - 1;
194            for (int i = 0; i < nd2; i++) {
195                tmp = a[offa + i];
196                a[offa + i] = a[idx0 - i];
197                a[idx0 - i] = tmp;
198            }
199        }
200        dct.inverse(a, offa, scale);
201        int startidx = 1 + offa;
202        int stopidx = offa + n;
203        for (int i = startidx; i < stopidx; i += 2) {
204            a[i] = -a[i];
205        }
206    }
207}