comparison libs/commons-math-2.1/docs/apidocs/src-html/org/apache/commons/math/transform/FastFourierTransformer.html @ 13:cbf34dd4d7e6

commons-math-2.1 added
author dwinter
date Tue, 04 Jan 2011 10:02:07 +0100
parents
children
comparison
equal deleted inserted replaced
12:970d26a94fb7 13:cbf34dd4d7e6
1 <HTML>
2 <BODY BGCOLOR="white">
3 <PRE>
4 <FONT color="green">001</FONT> /*<a name="line.1"></a>
5 <FONT color="green">002</FONT> * Licensed to the Apache Software Foundation (ASF) under one or more<a name="line.2"></a>
6 <FONT color="green">003</FONT> * contributor license agreements. See the NOTICE file distributed with<a name="line.3"></a>
7 <FONT color="green">004</FONT> * this work for additional information regarding copyright ownership.<a name="line.4"></a>
8 <FONT color="green">005</FONT> * The ASF licenses this file to You under the Apache License, Version 2.0<a name="line.5"></a>
9 <FONT color="green">006</FONT> * (the "License"); you may not use this file except in compliance with<a name="line.6"></a>
10 <FONT color="green">007</FONT> * the License. You may obtain a copy of the License at<a name="line.7"></a>
11 <FONT color="green">008</FONT> *<a name="line.8"></a>
12 <FONT color="green">009</FONT> * http://www.apache.org/licenses/LICENSE-2.0<a name="line.9"></a>
13 <FONT color="green">010</FONT> *<a name="line.10"></a>
14 <FONT color="green">011</FONT> * Unless required by applicable law or agreed to in writing, software<a name="line.11"></a>
15 <FONT color="green">012</FONT> * distributed under the License is distributed on an "AS IS" BASIS,<a name="line.12"></a>
16 <FONT color="green">013</FONT> * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.<a name="line.13"></a>
17 <FONT color="green">014</FONT> * See the License for the specific language governing permissions and<a name="line.14"></a>
18 <FONT color="green">015</FONT> * limitations under the License.<a name="line.15"></a>
19 <FONT color="green">016</FONT> */<a name="line.16"></a>
20 <FONT color="green">017</FONT> package org.apache.commons.math.transform;<a name="line.17"></a>
21 <FONT color="green">018</FONT> <a name="line.18"></a>
22 <FONT color="green">019</FONT> import java.io.Serializable;<a name="line.19"></a>
23 <FONT color="green">020</FONT> import java.lang.reflect.Array;<a name="line.20"></a>
24 <FONT color="green">021</FONT> <a name="line.21"></a>
25 <FONT color="green">022</FONT> import org.apache.commons.math.FunctionEvaluationException;<a name="line.22"></a>
26 <FONT color="green">023</FONT> import org.apache.commons.math.MathRuntimeException;<a name="line.23"></a>
27 <FONT color="green">024</FONT> import org.apache.commons.math.analysis.UnivariateRealFunction;<a name="line.24"></a>
28 <FONT color="green">025</FONT> import org.apache.commons.math.complex.Complex;<a name="line.25"></a>
29 <FONT color="green">026</FONT> <a name="line.26"></a>
30 <FONT color="green">027</FONT> /**<a name="line.27"></a>
31 <FONT color="green">028</FONT> * Implements the &lt;a href="http://mathworld.wolfram.com/FastFourierTransform.html"&gt;<a name="line.28"></a>
32 <FONT color="green">029</FONT> * Fast Fourier Transform&lt;/a&gt; for transformation of one-dimensional data sets.<a name="line.29"></a>
33 <FONT color="green">030</FONT> * For reference, see &lt;b&gt;Applied Numerical Linear Algebra&lt;/b&gt;, ISBN 0898713897,<a name="line.30"></a>
34 <FONT color="green">031</FONT> * chapter 6.<a name="line.31"></a>
35 <FONT color="green">032</FONT> * &lt;p&gt;<a name="line.32"></a>
36 <FONT color="green">033</FONT> * There are several conventions for the definition of FFT and inverse FFT,<a name="line.33"></a>
37 <FONT color="green">034</FONT> * mainly on different coefficient and exponent. Here the equations are listed<a name="line.34"></a>
38 <FONT color="green">035</FONT> * in the comments of the corresponding methods.&lt;/p&gt;<a name="line.35"></a>
39 <FONT color="green">036</FONT> * &lt;p&gt;<a name="line.36"></a>
40 <FONT color="green">037</FONT> * We require the length of data set to be power of 2, this greatly simplifies<a name="line.37"></a>
41 <FONT color="green">038</FONT> * and speeds up the code. Users can pad the data with zeros to meet this<a name="line.38"></a>
42 <FONT color="green">039</FONT> * requirement. There are other flavors of FFT, for reference, see S. Winograd,<a name="line.39"></a>
43 <FONT color="green">040</FONT> * &lt;i&gt;On computing the discrete Fourier transform&lt;/i&gt;, Mathematics of Computation,<a name="line.40"></a>
44 <FONT color="green">041</FONT> * 32 (1978), 175 - 199.&lt;/p&gt;<a name="line.41"></a>
45 <FONT color="green">042</FONT> *<a name="line.42"></a>
46 <FONT color="green">043</FONT> * @version $Revision: 885278 $ $Date: 2009-11-29 16:47:51 -0500 (Sun, 29 Nov 2009) $<a name="line.43"></a>
47 <FONT color="green">044</FONT> * @since 1.2<a name="line.44"></a>
48 <FONT color="green">045</FONT> */<a name="line.45"></a>
49 <FONT color="green">046</FONT> public class FastFourierTransformer implements Serializable {<a name="line.46"></a>
50 <FONT color="green">047</FONT> <a name="line.47"></a>
51 <FONT color="green">048</FONT> /** Serializable version identifier. */<a name="line.48"></a>
52 <FONT color="green">049</FONT> static final long serialVersionUID = 5138259215438106000L;<a name="line.49"></a>
53 <FONT color="green">050</FONT> <a name="line.50"></a>
54 <FONT color="green">051</FONT> /** Message for not power of 2. */<a name="line.51"></a>
55 <FONT color="green">052</FONT> private static final String NOT_POWER_OF_TWO_MESSAGE =<a name="line.52"></a>
56 <FONT color="green">053</FONT> "{0} is not a power of 2, consider padding for fix";<a name="line.53"></a>
57 <FONT color="green">054</FONT> <a name="line.54"></a>
58 <FONT color="green">055</FONT> /** Message for dimension mismatch. */<a name="line.55"></a>
59 <FONT color="green">056</FONT> private static final String DIMENSION_MISMATCH_MESSAGE =<a name="line.56"></a>
60 <FONT color="green">057</FONT> "some dimensions don't match: {0} != {1}";<a name="line.57"></a>
61 <FONT color="green">058</FONT> <a name="line.58"></a>
62 <FONT color="green">059</FONT> /** Message for not computed roots of unity. */<a name="line.59"></a>
63 <FONT color="green">060</FONT> private static final String MISSING_ROOTS_OF_UNITY_MESSAGE =<a name="line.60"></a>
64 <FONT color="green">061</FONT> "roots of unity have not been computed yet";<a name="line.61"></a>
65 <FONT color="green">062</FONT> <a name="line.62"></a>
66 <FONT color="green">063</FONT> /** Message for out of range root index. */<a name="line.63"></a>
67 <FONT color="green">064</FONT> private static final String OUT_OF_RANGE_ROOT_INDEX_MESSAGE =<a name="line.64"></a>
68 <FONT color="green">065</FONT> "out of range root of unity index {0} (must be in [{1};{2}])";<a name="line.65"></a>
69 <FONT color="green">066</FONT> <a name="line.66"></a>
70 <FONT color="green">067</FONT> /** roots of unity */<a name="line.67"></a>
71 <FONT color="green">068</FONT> private RootsOfUnity roots = new RootsOfUnity();<a name="line.68"></a>
72 <FONT color="green">069</FONT> <a name="line.69"></a>
73 <FONT color="green">070</FONT> /**<a name="line.70"></a>
74 <FONT color="green">071</FONT> * Construct a default transformer.<a name="line.71"></a>
75 <FONT color="green">072</FONT> */<a name="line.72"></a>
76 <FONT color="green">073</FONT> public FastFourierTransformer() {<a name="line.73"></a>
77 <FONT color="green">074</FONT> super();<a name="line.74"></a>
78 <FONT color="green">075</FONT> }<a name="line.75"></a>
79 <FONT color="green">076</FONT> <a name="line.76"></a>
80 <FONT color="green">077</FONT> /**<a name="line.77"></a>
81 <FONT color="green">078</FONT> * Transform the given real data set.<a name="line.78"></a>
82 <FONT color="green">079</FONT> * &lt;p&gt;<a name="line.79"></a>
83 <FONT color="green">080</FONT> * The formula is $ y_n = \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k $<a name="line.80"></a>
84 <FONT color="green">081</FONT> * &lt;/p&gt;<a name="line.81"></a>
85 <FONT color="green">082</FONT> *<a name="line.82"></a>
86 <FONT color="green">083</FONT> * @param f the real data array to be transformed<a name="line.83"></a>
87 <FONT color="green">084</FONT> * @return the complex transformed array<a name="line.84"></a>
88 <FONT color="green">085</FONT> * @throws IllegalArgumentException if any parameters are invalid<a name="line.85"></a>
89 <FONT color="green">086</FONT> */<a name="line.86"></a>
90 <FONT color="green">087</FONT> public Complex[] transform(double f[])<a name="line.87"></a>
91 <FONT color="green">088</FONT> throws IllegalArgumentException {<a name="line.88"></a>
92 <FONT color="green">089</FONT> return fft(f, false);<a name="line.89"></a>
93 <FONT color="green">090</FONT> }<a name="line.90"></a>
94 <FONT color="green">091</FONT> <a name="line.91"></a>
95 <FONT color="green">092</FONT> /**<a name="line.92"></a>
96 <FONT color="green">093</FONT> * Transform the given real function, sampled on the given interval.<a name="line.93"></a>
97 <FONT color="green">094</FONT> * &lt;p&gt;<a name="line.94"></a>
98 <FONT color="green">095</FONT> * The formula is $ y_n = \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k $<a name="line.95"></a>
99 <FONT color="green">096</FONT> * &lt;/p&gt;<a name="line.96"></a>
100 <FONT color="green">097</FONT> *<a name="line.97"></a>
101 <FONT color="green">098</FONT> * @param f the function to be sampled and transformed<a name="line.98"></a>
102 <FONT color="green">099</FONT> * @param min the lower bound for the interval<a name="line.99"></a>
103 <FONT color="green">100</FONT> * @param max the upper bound for the interval<a name="line.100"></a>
104 <FONT color="green">101</FONT> * @param n the number of sample points<a name="line.101"></a>
105 <FONT color="green">102</FONT> * @return the complex transformed array<a name="line.102"></a>
106 <FONT color="green">103</FONT> * @throws FunctionEvaluationException if function cannot be evaluated<a name="line.103"></a>
107 <FONT color="green">104</FONT> * at some point<a name="line.104"></a>
108 <FONT color="green">105</FONT> * @throws IllegalArgumentException if any parameters are invalid<a name="line.105"></a>
109 <FONT color="green">106</FONT> */<a name="line.106"></a>
110 <FONT color="green">107</FONT> public Complex[] transform(UnivariateRealFunction f,<a name="line.107"></a>
111 <FONT color="green">108</FONT> double min, double max, int n)<a name="line.108"></a>
112 <FONT color="green">109</FONT> throws FunctionEvaluationException, IllegalArgumentException {<a name="line.109"></a>
113 <FONT color="green">110</FONT> double data[] = sample(f, min, max, n);<a name="line.110"></a>
114 <FONT color="green">111</FONT> return fft(data, false);<a name="line.111"></a>
115 <FONT color="green">112</FONT> }<a name="line.112"></a>
116 <FONT color="green">113</FONT> <a name="line.113"></a>
117 <FONT color="green">114</FONT> /**<a name="line.114"></a>
118 <FONT color="green">115</FONT> * Transform the given complex data set.<a name="line.115"></a>
119 <FONT color="green">116</FONT> * &lt;p&gt;<a name="line.116"></a>
120 <FONT color="green">117</FONT> * The formula is $ y_n = \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k $<a name="line.117"></a>
121 <FONT color="green">118</FONT> * &lt;/p&gt;<a name="line.118"></a>
122 <FONT color="green">119</FONT> *<a name="line.119"></a>
123 <FONT color="green">120</FONT> * @param f the complex data array to be transformed<a name="line.120"></a>
124 <FONT color="green">121</FONT> * @return the complex transformed array<a name="line.121"></a>
125 <FONT color="green">122</FONT> * @throws IllegalArgumentException if any parameters are invalid<a name="line.122"></a>
126 <FONT color="green">123</FONT> */<a name="line.123"></a>
127 <FONT color="green">124</FONT> public Complex[] transform(Complex f[])<a name="line.124"></a>
128 <FONT color="green">125</FONT> throws IllegalArgumentException {<a name="line.125"></a>
129 <FONT color="green">126</FONT> roots.computeOmega(f.length);<a name="line.126"></a>
130 <FONT color="green">127</FONT> return fft(f);<a name="line.127"></a>
131 <FONT color="green">128</FONT> }<a name="line.128"></a>
132 <FONT color="green">129</FONT> <a name="line.129"></a>
133 <FONT color="green">130</FONT> /**<a name="line.130"></a>
134 <FONT color="green">131</FONT> * Transform the given real data set.<a name="line.131"></a>
135 <FONT color="green">132</FONT> * &lt;p&gt;<a name="line.132"></a>
136 <FONT color="green">133</FONT> * The formula is $y_n = (1/\sqrt{N}) \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k$<a name="line.133"></a>
137 <FONT color="green">134</FONT> * &lt;/p&gt;<a name="line.134"></a>
138 <FONT color="green">135</FONT> *<a name="line.135"></a>
139 <FONT color="green">136</FONT> * @param f the real data array to be transformed<a name="line.136"></a>
140 <FONT color="green">137</FONT> * @return the complex transformed array<a name="line.137"></a>
141 <FONT color="green">138</FONT> * @throws IllegalArgumentException if any parameters are invalid<a name="line.138"></a>
142 <FONT color="green">139</FONT> */<a name="line.139"></a>
143 <FONT color="green">140</FONT> public Complex[] transform2(double f[])<a name="line.140"></a>
144 <FONT color="green">141</FONT> throws IllegalArgumentException {<a name="line.141"></a>
145 <FONT color="green">142</FONT> <a name="line.142"></a>
146 <FONT color="green">143</FONT> double scaling_coefficient = 1.0 / Math.sqrt(f.length);<a name="line.143"></a>
147 <FONT color="green">144</FONT> return scaleArray(fft(f, false), scaling_coefficient);<a name="line.144"></a>
148 <FONT color="green">145</FONT> }<a name="line.145"></a>
149 <FONT color="green">146</FONT> <a name="line.146"></a>
150 <FONT color="green">147</FONT> /**<a name="line.147"></a>
151 <FONT color="green">148</FONT> * Transform the given real function, sampled on the given interval.<a name="line.148"></a>
152 <FONT color="green">149</FONT> * &lt;p&gt;<a name="line.149"></a>
153 <FONT color="green">150</FONT> * The formula is $y_n = (1/\sqrt{N}) \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k$<a name="line.150"></a>
154 <FONT color="green">151</FONT> * &lt;/p&gt;<a name="line.151"></a>
155 <FONT color="green">152</FONT> *<a name="line.152"></a>
156 <FONT color="green">153</FONT> * @param f the function to be sampled and transformed<a name="line.153"></a>
157 <FONT color="green">154</FONT> * @param min the lower bound for the interval<a name="line.154"></a>
158 <FONT color="green">155</FONT> * @param max the upper bound for the interval<a name="line.155"></a>
159 <FONT color="green">156</FONT> * @param n the number of sample points<a name="line.156"></a>
160 <FONT color="green">157</FONT> * @return the complex transformed array<a name="line.157"></a>
161 <FONT color="green">158</FONT> * @throws FunctionEvaluationException if function cannot be evaluated<a name="line.158"></a>
162 <FONT color="green">159</FONT> * at some point<a name="line.159"></a>
163 <FONT color="green">160</FONT> * @throws IllegalArgumentException if any parameters are invalid<a name="line.160"></a>
164 <FONT color="green">161</FONT> */<a name="line.161"></a>
165 <FONT color="green">162</FONT> public Complex[] transform2(UnivariateRealFunction f,<a name="line.162"></a>
166 <FONT color="green">163</FONT> double min, double max, int n)<a name="line.163"></a>
167 <FONT color="green">164</FONT> throws FunctionEvaluationException, IllegalArgumentException {<a name="line.164"></a>
168 <FONT color="green">165</FONT> <a name="line.165"></a>
169 <FONT color="green">166</FONT> double data[] = sample(f, min, max, n);<a name="line.166"></a>
170 <FONT color="green">167</FONT> double scaling_coefficient = 1.0 / Math.sqrt(n);<a name="line.167"></a>
171 <FONT color="green">168</FONT> return scaleArray(fft(data, false), scaling_coefficient);<a name="line.168"></a>
172 <FONT color="green">169</FONT> }<a name="line.169"></a>
173 <FONT color="green">170</FONT> <a name="line.170"></a>
174 <FONT color="green">171</FONT> /**<a name="line.171"></a>
175 <FONT color="green">172</FONT> * Transform the given complex data set.<a name="line.172"></a>
176 <FONT color="green">173</FONT> * &lt;p&gt;<a name="line.173"></a>
177 <FONT color="green">174</FONT> * The formula is $y_n = (1/\sqrt{N}) \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k$<a name="line.174"></a>
178 <FONT color="green">175</FONT> * &lt;/p&gt;<a name="line.175"></a>
179 <FONT color="green">176</FONT> *<a name="line.176"></a>
180 <FONT color="green">177</FONT> * @param f the complex data array to be transformed<a name="line.177"></a>
181 <FONT color="green">178</FONT> * @return the complex transformed array<a name="line.178"></a>
182 <FONT color="green">179</FONT> * @throws IllegalArgumentException if any parameters are invalid<a name="line.179"></a>
183 <FONT color="green">180</FONT> */<a name="line.180"></a>
184 <FONT color="green">181</FONT> public Complex[] transform2(Complex f[])<a name="line.181"></a>
185 <FONT color="green">182</FONT> throws IllegalArgumentException {<a name="line.182"></a>
186 <FONT color="green">183</FONT> <a name="line.183"></a>
187 <FONT color="green">184</FONT> roots.computeOmega(f.length);<a name="line.184"></a>
188 <FONT color="green">185</FONT> double scaling_coefficient = 1.0 / Math.sqrt(f.length);<a name="line.185"></a>
189 <FONT color="green">186</FONT> return scaleArray(fft(f), scaling_coefficient);<a name="line.186"></a>
190 <FONT color="green">187</FONT> }<a name="line.187"></a>
191 <FONT color="green">188</FONT> <a name="line.188"></a>
192 <FONT color="green">189</FONT> /**<a name="line.189"></a>
193 <FONT color="green">190</FONT> * Inversely transform the given real data set.<a name="line.190"></a>
194 <FONT color="green">191</FONT> * &lt;p&gt;<a name="line.191"></a>
195 <FONT color="green">192</FONT> * The formula is $ x_k = (1/N) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n $<a name="line.192"></a>
196 <FONT color="green">193</FONT> * &lt;/p&gt;<a name="line.193"></a>
197 <FONT color="green">194</FONT> *<a name="line.194"></a>
198 <FONT color="green">195</FONT> * @param f the real data array to be inversely transformed<a name="line.195"></a>
199 <FONT color="green">196</FONT> * @return the complex inversely transformed array<a name="line.196"></a>
200 <FONT color="green">197</FONT> * @throws IllegalArgumentException if any parameters are invalid<a name="line.197"></a>
201 <FONT color="green">198</FONT> */<a name="line.198"></a>
202 <FONT color="green">199</FONT> public Complex[] inversetransform(double f[])<a name="line.199"></a>
203 <FONT color="green">200</FONT> throws IllegalArgumentException {<a name="line.200"></a>
204 <FONT color="green">201</FONT> <a name="line.201"></a>
205 <FONT color="green">202</FONT> double scaling_coefficient = 1.0 / f.length;<a name="line.202"></a>
206 <FONT color="green">203</FONT> return scaleArray(fft(f, true), scaling_coefficient);<a name="line.203"></a>
207 <FONT color="green">204</FONT> }<a name="line.204"></a>
208 <FONT color="green">205</FONT> <a name="line.205"></a>
209 <FONT color="green">206</FONT> /**<a name="line.206"></a>
210 <FONT color="green">207</FONT> * Inversely transform the given real function, sampled on the given interval.<a name="line.207"></a>
211 <FONT color="green">208</FONT> * &lt;p&gt;<a name="line.208"></a>
212 <FONT color="green">209</FONT> * The formula is $ x_k = (1/N) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n $<a name="line.209"></a>
213 <FONT color="green">210</FONT> * &lt;/p&gt;<a name="line.210"></a>
214 <FONT color="green">211</FONT> *<a name="line.211"></a>
215 <FONT color="green">212</FONT> * @param f the function to be sampled and inversely transformed<a name="line.212"></a>
216 <FONT color="green">213</FONT> * @param min the lower bound for the interval<a name="line.213"></a>
217 <FONT color="green">214</FONT> * @param max the upper bound for the interval<a name="line.214"></a>
218 <FONT color="green">215</FONT> * @param n the number of sample points<a name="line.215"></a>
219 <FONT color="green">216</FONT> * @return the complex inversely transformed array<a name="line.216"></a>
220 <FONT color="green">217</FONT> * @throws FunctionEvaluationException if function cannot be evaluated<a name="line.217"></a>
221 <FONT color="green">218</FONT> * at some point<a name="line.218"></a>
222 <FONT color="green">219</FONT> * @throws IllegalArgumentException if any parameters are invalid<a name="line.219"></a>
223 <FONT color="green">220</FONT> */<a name="line.220"></a>
224 <FONT color="green">221</FONT> public Complex[] inversetransform(UnivariateRealFunction f,<a name="line.221"></a>
225 <FONT color="green">222</FONT> double min, double max, int n)<a name="line.222"></a>
226 <FONT color="green">223</FONT> throws FunctionEvaluationException, IllegalArgumentException {<a name="line.223"></a>
227 <FONT color="green">224</FONT> <a name="line.224"></a>
228 <FONT color="green">225</FONT> double data[] = sample(f, min, max, n);<a name="line.225"></a>
229 <FONT color="green">226</FONT> double scaling_coefficient = 1.0 / n;<a name="line.226"></a>
230 <FONT color="green">227</FONT> return scaleArray(fft(data, true), scaling_coefficient);<a name="line.227"></a>
231 <FONT color="green">228</FONT> }<a name="line.228"></a>
232 <FONT color="green">229</FONT> <a name="line.229"></a>
233 <FONT color="green">230</FONT> /**<a name="line.230"></a>
234 <FONT color="green">231</FONT> * Inversely transform the given complex data set.<a name="line.231"></a>
235 <FONT color="green">232</FONT> * &lt;p&gt;<a name="line.232"></a>
236 <FONT color="green">233</FONT> * The formula is $ x_k = (1/N) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n $<a name="line.233"></a>
237 <FONT color="green">234</FONT> * &lt;/p&gt;<a name="line.234"></a>
238 <FONT color="green">235</FONT> *<a name="line.235"></a>
239 <FONT color="green">236</FONT> * @param f the complex data array to be inversely transformed<a name="line.236"></a>
240 <FONT color="green">237</FONT> * @return the complex inversely transformed array<a name="line.237"></a>
241 <FONT color="green">238</FONT> * @throws IllegalArgumentException if any parameters are invalid<a name="line.238"></a>
242 <FONT color="green">239</FONT> */<a name="line.239"></a>
243 <FONT color="green">240</FONT> public Complex[] inversetransform(Complex f[])<a name="line.240"></a>
244 <FONT color="green">241</FONT> throws IllegalArgumentException {<a name="line.241"></a>
245 <FONT color="green">242</FONT> <a name="line.242"></a>
246 <FONT color="green">243</FONT> roots.computeOmega(-f.length); // pass negative argument<a name="line.243"></a>
247 <FONT color="green">244</FONT> double scaling_coefficient = 1.0 / f.length;<a name="line.244"></a>
248 <FONT color="green">245</FONT> return scaleArray(fft(f), scaling_coefficient);<a name="line.245"></a>
249 <FONT color="green">246</FONT> }<a name="line.246"></a>
250 <FONT color="green">247</FONT> <a name="line.247"></a>
251 <FONT color="green">248</FONT> /**<a name="line.248"></a>
252 <FONT color="green">249</FONT> * Inversely transform the given real data set.<a name="line.249"></a>
253 <FONT color="green">250</FONT> * &lt;p&gt;<a name="line.250"></a>
254 <FONT color="green">251</FONT> * The formula is $x_k = (1/\sqrt{N}) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n$<a name="line.251"></a>
255 <FONT color="green">252</FONT> * &lt;/p&gt;<a name="line.252"></a>
256 <FONT color="green">253</FONT> *<a name="line.253"></a>
257 <FONT color="green">254</FONT> * @param f the real data array to be inversely transformed<a name="line.254"></a>
258 <FONT color="green">255</FONT> * @return the complex inversely transformed array<a name="line.255"></a>
259 <FONT color="green">256</FONT> * @throws IllegalArgumentException if any parameters are invalid<a name="line.256"></a>
260 <FONT color="green">257</FONT> */<a name="line.257"></a>
261 <FONT color="green">258</FONT> public Complex[] inversetransform2(double f[])<a name="line.258"></a>
262 <FONT color="green">259</FONT> throws IllegalArgumentException {<a name="line.259"></a>
263 <FONT color="green">260</FONT> <a name="line.260"></a>
264 <FONT color="green">261</FONT> double scaling_coefficient = 1.0 / Math.sqrt(f.length);<a name="line.261"></a>
265 <FONT color="green">262</FONT> return scaleArray(fft(f, true), scaling_coefficient);<a name="line.262"></a>
266 <FONT color="green">263</FONT> }<a name="line.263"></a>
267 <FONT color="green">264</FONT> <a name="line.264"></a>
268 <FONT color="green">265</FONT> /**<a name="line.265"></a>
269 <FONT color="green">266</FONT> * Inversely transform the given real function, sampled on the given interval.<a name="line.266"></a>
270 <FONT color="green">267</FONT> * &lt;p&gt;<a name="line.267"></a>
271 <FONT color="green">268</FONT> * The formula is $x_k = (1/\sqrt{N}) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n$<a name="line.268"></a>
272 <FONT color="green">269</FONT> * &lt;/p&gt;<a name="line.269"></a>
273 <FONT color="green">270</FONT> *<a name="line.270"></a>
274 <FONT color="green">271</FONT> * @param f the function to be sampled and inversely transformed<a name="line.271"></a>
275 <FONT color="green">272</FONT> * @param min the lower bound for the interval<a name="line.272"></a>
276 <FONT color="green">273</FONT> * @param max the upper bound for the interval<a name="line.273"></a>
277 <FONT color="green">274</FONT> * @param n the number of sample points<a name="line.274"></a>
278 <FONT color="green">275</FONT> * @return the complex inversely transformed array<a name="line.275"></a>
279 <FONT color="green">276</FONT> * @throws FunctionEvaluationException if function cannot be evaluated<a name="line.276"></a>
280 <FONT color="green">277</FONT> * at some point<a name="line.277"></a>
281 <FONT color="green">278</FONT> * @throws IllegalArgumentException if any parameters are invalid<a name="line.278"></a>
282 <FONT color="green">279</FONT> */<a name="line.279"></a>
283 <FONT color="green">280</FONT> public Complex[] inversetransform2(UnivariateRealFunction f,<a name="line.280"></a>
284 <FONT color="green">281</FONT> double min, double max, int n)<a name="line.281"></a>
285 <FONT color="green">282</FONT> throws FunctionEvaluationException, IllegalArgumentException {<a name="line.282"></a>
286 <FONT color="green">283</FONT> <a name="line.283"></a>
287 <FONT color="green">284</FONT> double data[] = sample(f, min, max, n);<a name="line.284"></a>
288 <FONT color="green">285</FONT> double scaling_coefficient = 1.0 / Math.sqrt(n);<a name="line.285"></a>
289 <FONT color="green">286</FONT> return scaleArray(fft(data, true), scaling_coefficient);<a name="line.286"></a>
290 <FONT color="green">287</FONT> }<a name="line.287"></a>
291 <FONT color="green">288</FONT> <a name="line.288"></a>
292 <FONT color="green">289</FONT> /**<a name="line.289"></a>
293 <FONT color="green">290</FONT> * Inversely transform the given complex data set.<a name="line.290"></a>
294 <FONT color="green">291</FONT> * &lt;p&gt;<a name="line.291"></a>
295 <FONT color="green">292</FONT> * The formula is $x_k = (1/\sqrt{N}) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n$<a name="line.292"></a>
296 <FONT color="green">293</FONT> * &lt;/p&gt;<a name="line.293"></a>
297 <FONT color="green">294</FONT> *<a name="line.294"></a>
298 <FONT color="green">295</FONT> * @param f the complex data array to be inversely transformed<a name="line.295"></a>
299 <FONT color="green">296</FONT> * @return the complex inversely transformed array<a name="line.296"></a>
300 <FONT color="green">297</FONT> * @throws IllegalArgumentException if any parameters are invalid<a name="line.297"></a>
301 <FONT color="green">298</FONT> */<a name="line.298"></a>
302 <FONT color="green">299</FONT> public Complex[] inversetransform2(Complex f[])<a name="line.299"></a>
303 <FONT color="green">300</FONT> throws IllegalArgumentException {<a name="line.300"></a>
304 <FONT color="green">301</FONT> <a name="line.301"></a>
305 <FONT color="green">302</FONT> roots.computeOmega(-f.length); // pass negative argument<a name="line.302"></a>
306 <FONT color="green">303</FONT> double scaling_coefficient = 1.0 / Math.sqrt(f.length);<a name="line.303"></a>
307 <FONT color="green">304</FONT> return scaleArray(fft(f), scaling_coefficient);<a name="line.304"></a>
308 <FONT color="green">305</FONT> }<a name="line.305"></a>
309 <FONT color="green">306</FONT> <a name="line.306"></a>
310 <FONT color="green">307</FONT> /**<a name="line.307"></a>
311 <FONT color="green">308</FONT> * Perform the base-4 Cooley-Tukey FFT algorithm (including inverse).<a name="line.308"></a>
312 <FONT color="green">309</FONT> *<a name="line.309"></a>
313 <FONT color="green">310</FONT> * @param f the real data array to be transformed<a name="line.310"></a>
314 <FONT color="green">311</FONT> * @param isInverse the indicator of forward or inverse transform<a name="line.311"></a>
315 <FONT color="green">312</FONT> * @return the complex transformed array<a name="line.312"></a>
316 <FONT color="green">313</FONT> * @throws IllegalArgumentException if any parameters are invalid<a name="line.313"></a>
317 <FONT color="green">314</FONT> */<a name="line.314"></a>
318 <FONT color="green">315</FONT> protected Complex[] fft(double f[], boolean isInverse)<a name="line.315"></a>
319 <FONT color="green">316</FONT> throws IllegalArgumentException {<a name="line.316"></a>
320 <FONT color="green">317</FONT> <a name="line.317"></a>
321 <FONT color="green">318</FONT> verifyDataSet(f);<a name="line.318"></a>
322 <FONT color="green">319</FONT> Complex F[] = new Complex[f.length];<a name="line.319"></a>
323 <FONT color="green">320</FONT> if (f.length == 1) {<a name="line.320"></a>
324 <FONT color="green">321</FONT> F[0] = new Complex(f[0], 0.0);<a name="line.321"></a>
325 <FONT color="green">322</FONT> return F;<a name="line.322"></a>
326 <FONT color="green">323</FONT> }<a name="line.323"></a>
327 <FONT color="green">324</FONT> <a name="line.324"></a>
328 <FONT color="green">325</FONT> // Rather than the naive real to complex conversion, pack 2N<a name="line.325"></a>
329 <FONT color="green">326</FONT> // real numbers into N complex numbers for better performance.<a name="line.326"></a>
330 <FONT color="green">327</FONT> int N = f.length &gt;&gt; 1;<a name="line.327"></a>
331 <FONT color="green">328</FONT> Complex c[] = new Complex[N];<a name="line.328"></a>
332 <FONT color="green">329</FONT> for (int i = 0; i &lt; N; i++) {<a name="line.329"></a>
333 <FONT color="green">330</FONT> c[i] = new Complex(f[2*i], f[2*i+1]);<a name="line.330"></a>
334 <FONT color="green">331</FONT> }<a name="line.331"></a>
335 <FONT color="green">332</FONT> roots.computeOmega(isInverse ? -N : N);<a name="line.332"></a>
336 <FONT color="green">333</FONT> Complex z[] = fft(c);<a name="line.333"></a>
337 <FONT color="green">334</FONT> <a name="line.334"></a>
338 <FONT color="green">335</FONT> // reconstruct the FFT result for the original array<a name="line.335"></a>
339 <FONT color="green">336</FONT> roots.computeOmega(isInverse ? -2*N : 2*N);<a name="line.336"></a>
340 <FONT color="green">337</FONT> F[0] = new Complex(2 * (z[0].getReal() + z[0].getImaginary()), 0.0);<a name="line.337"></a>
341 <FONT color="green">338</FONT> F[N] = new Complex(2 * (z[0].getReal() - z[0].getImaginary()), 0.0);<a name="line.338"></a>
342 <FONT color="green">339</FONT> for (int i = 1; i &lt; N; i++) {<a name="line.339"></a>
343 <FONT color="green">340</FONT> Complex A = z[N-i].conjugate();<a name="line.340"></a>
344 <FONT color="green">341</FONT> Complex B = z[i].add(A);<a name="line.341"></a>
345 <FONT color="green">342</FONT> Complex C = z[i].subtract(A);<a name="line.342"></a>
346 <FONT color="green">343</FONT> //Complex D = roots.getOmega(i).multiply(Complex.I);<a name="line.343"></a>
347 <FONT color="green">344</FONT> Complex D = new Complex(-roots.getOmegaImaginary(i),<a name="line.344"></a>
348 <FONT color="green">345</FONT> roots.getOmegaReal(i));<a name="line.345"></a>
349 <FONT color="green">346</FONT> F[i] = B.subtract(C.multiply(D));<a name="line.346"></a>
350 <FONT color="green">347</FONT> F[2*N-i] = F[i].conjugate();<a name="line.347"></a>
351 <FONT color="green">348</FONT> }<a name="line.348"></a>
352 <FONT color="green">349</FONT> <a name="line.349"></a>
353 <FONT color="green">350</FONT> return scaleArray(F, 0.5);<a name="line.350"></a>
354 <FONT color="green">351</FONT> }<a name="line.351"></a>
355 <FONT color="green">352</FONT> <a name="line.352"></a>
356 <FONT color="green">353</FONT> /**<a name="line.353"></a>
357 <FONT color="green">354</FONT> * Perform the base-4 Cooley-Tukey FFT algorithm (including inverse).<a name="line.354"></a>
358 <FONT color="green">355</FONT> *<a name="line.355"></a>
359 <FONT color="green">356</FONT> * @param data the complex data array to be transformed<a name="line.356"></a>
360 <FONT color="green">357</FONT> * @return the complex transformed array<a name="line.357"></a>
361 <FONT color="green">358</FONT> * @throws IllegalArgumentException if any parameters are invalid<a name="line.358"></a>
362 <FONT color="green">359</FONT> */<a name="line.359"></a>
363 <FONT color="green">360</FONT> protected Complex[] fft(Complex data[])<a name="line.360"></a>
364 <FONT color="green">361</FONT> throws IllegalArgumentException {<a name="line.361"></a>
365 <FONT color="green">362</FONT> <a name="line.362"></a>
366 <FONT color="green">363</FONT> final int n = data.length;<a name="line.363"></a>
367 <FONT color="green">364</FONT> final Complex f[] = new Complex[n];<a name="line.364"></a>
368 <FONT color="green">365</FONT> <a name="line.365"></a>
369 <FONT color="green">366</FONT> // initial simple cases<a name="line.366"></a>
370 <FONT color="green">367</FONT> verifyDataSet(data);<a name="line.367"></a>
371 <FONT color="green">368</FONT> if (n == 1) {<a name="line.368"></a>
372 <FONT color="green">369</FONT> f[0] = data[0];<a name="line.369"></a>
373 <FONT color="green">370</FONT> return f;<a name="line.370"></a>
374 <FONT color="green">371</FONT> }<a name="line.371"></a>
375 <FONT color="green">372</FONT> if (n == 2) {<a name="line.372"></a>
376 <FONT color="green">373</FONT> f[0] = data[0].add(data[1]);<a name="line.373"></a>
377 <FONT color="green">374</FONT> f[1] = data[0].subtract(data[1]);<a name="line.374"></a>
378 <FONT color="green">375</FONT> return f;<a name="line.375"></a>
379 <FONT color="green">376</FONT> }<a name="line.376"></a>
380 <FONT color="green">377</FONT> <a name="line.377"></a>
381 <FONT color="green">378</FONT> // permute original data array in bit-reversal order<a name="line.378"></a>
382 <FONT color="green">379</FONT> int ii = 0;<a name="line.379"></a>
383 <FONT color="green">380</FONT> for (int i = 0; i &lt; n; i++) {<a name="line.380"></a>
384 <FONT color="green">381</FONT> f[i] = data[ii];<a name="line.381"></a>
385 <FONT color="green">382</FONT> int k = n &gt;&gt; 1;<a name="line.382"></a>
386 <FONT color="green">383</FONT> while (ii &gt;= k &amp;&amp; k &gt; 0) {<a name="line.383"></a>
387 <FONT color="green">384</FONT> ii -= k; k &gt;&gt;= 1;<a name="line.384"></a>
388 <FONT color="green">385</FONT> }<a name="line.385"></a>
389 <FONT color="green">386</FONT> ii += k;<a name="line.386"></a>
390 <FONT color="green">387</FONT> }<a name="line.387"></a>
391 <FONT color="green">388</FONT> <a name="line.388"></a>
392 <FONT color="green">389</FONT> // the bottom base-4 round<a name="line.389"></a>
393 <FONT color="green">390</FONT> for (int i = 0; i &lt; n; i += 4) {<a name="line.390"></a>
394 <FONT color="green">391</FONT> final Complex a = f[i].add(f[i+1]);<a name="line.391"></a>
395 <FONT color="green">392</FONT> final Complex b = f[i+2].add(f[i+3]);<a name="line.392"></a>
396 <FONT color="green">393</FONT> final Complex c = f[i].subtract(f[i+1]);<a name="line.393"></a>
397 <FONT color="green">394</FONT> final Complex d = f[i+2].subtract(f[i+3]);<a name="line.394"></a>
398 <FONT color="green">395</FONT> final Complex e1 = c.add(d.multiply(Complex.I));<a name="line.395"></a>
399 <FONT color="green">396</FONT> final Complex e2 = c.subtract(d.multiply(Complex.I));<a name="line.396"></a>
400 <FONT color="green">397</FONT> f[i] = a.add(b);<a name="line.397"></a>
401 <FONT color="green">398</FONT> f[i+2] = a.subtract(b);<a name="line.398"></a>
402 <FONT color="green">399</FONT> // omegaCount indicates forward or inverse transform<a name="line.399"></a>
403 <FONT color="green">400</FONT> f[i+1] = roots.isForward() ? e2 : e1;<a name="line.400"></a>
404 <FONT color="green">401</FONT> f[i+3] = roots.isForward() ? e1 : e2;<a name="line.401"></a>
405 <FONT color="green">402</FONT> }<a name="line.402"></a>
406 <FONT color="green">403</FONT> <a name="line.403"></a>
407 <FONT color="green">404</FONT> // iterations from bottom to top take O(N*logN) time<a name="line.404"></a>
408 <FONT color="green">405</FONT> for (int i = 4; i &lt; n; i &lt;&lt;= 1) {<a name="line.405"></a>
409 <FONT color="green">406</FONT> final int m = n / (i&lt;&lt;1);<a name="line.406"></a>
410 <FONT color="green">407</FONT> for (int j = 0; j &lt; n; j += i&lt;&lt;1) {<a name="line.407"></a>
411 <FONT color="green">408</FONT> for (int k = 0; k &lt; i; k++) {<a name="line.408"></a>
412 <FONT color="green">409</FONT> //z = f[i+j+k].multiply(roots.getOmega(k*m));<a name="line.409"></a>
413 <FONT color="green">410</FONT> final int k_times_m = k*m;<a name="line.410"></a>
414 <FONT color="green">411</FONT> final double omega_k_times_m_real = roots.getOmegaReal(k_times_m);<a name="line.411"></a>
415 <FONT color="green">412</FONT> final double omega_k_times_m_imaginary = roots.getOmegaImaginary(k_times_m);<a name="line.412"></a>
416 <FONT color="green">413</FONT> //z = f[i+j+k].multiply(omega[k*m]);<a name="line.413"></a>
417 <FONT color="green">414</FONT> final Complex z = new Complex(<a name="line.414"></a>
418 <FONT color="green">415</FONT> f[i+j+k].getReal() * omega_k_times_m_real -<a name="line.415"></a>
419 <FONT color="green">416</FONT> f[i+j+k].getImaginary() * omega_k_times_m_imaginary,<a name="line.416"></a>
420 <FONT color="green">417</FONT> f[i+j+k].getReal() * omega_k_times_m_imaginary +<a name="line.417"></a>
421 <FONT color="green">418</FONT> f[i+j+k].getImaginary() * omega_k_times_m_real);<a name="line.418"></a>
422 <FONT color="green">419</FONT> <a name="line.419"></a>
423 <FONT color="green">420</FONT> f[i+j+k] = f[j+k].subtract(z);<a name="line.420"></a>
424 <FONT color="green">421</FONT> f[j+k] = f[j+k].add(z);<a name="line.421"></a>
425 <FONT color="green">422</FONT> }<a name="line.422"></a>
426 <FONT color="green">423</FONT> }<a name="line.423"></a>
427 <FONT color="green">424</FONT> }<a name="line.424"></a>
428 <FONT color="green">425</FONT> return f;<a name="line.425"></a>
429 <FONT color="green">426</FONT> }<a name="line.426"></a>
430 <FONT color="green">427</FONT> <a name="line.427"></a>
431 <FONT color="green">428</FONT> /**<a name="line.428"></a>
432 <FONT color="green">429</FONT> * Sample the given univariate real function on the given interval.<a name="line.429"></a>
433 <FONT color="green">430</FONT> * &lt;p&gt;<a name="line.430"></a>
434 <FONT color="green">431</FONT> * The interval is divided equally into N sections and sample points<a name="line.431"></a>
435 <FONT color="green">432</FONT> * are taken from min to max-(max-min)/N. Usually f(x) is periodic<a name="line.432"></a>
436 <FONT color="green">433</FONT> * such that f(min) = f(max) (note max is not sampled), but we don't<a name="line.433"></a>
437 <FONT color="green">434</FONT> * require that.&lt;/p&gt;<a name="line.434"></a>
438 <FONT color="green">435</FONT> *<a name="line.435"></a>
439 <FONT color="green">436</FONT> * @param f the function to be sampled<a name="line.436"></a>
440 <FONT color="green">437</FONT> * @param min the lower bound for the interval<a name="line.437"></a>
441 <FONT color="green">438</FONT> * @param max the upper bound for the interval<a name="line.438"></a>
442 <FONT color="green">439</FONT> * @param n the number of sample points<a name="line.439"></a>
443 <FONT color="green">440</FONT> * @return the samples array<a name="line.440"></a>
444 <FONT color="green">441</FONT> * @throws FunctionEvaluationException if function cannot be evaluated<a name="line.441"></a>
445 <FONT color="green">442</FONT> * at some point<a name="line.442"></a>
446 <FONT color="green">443</FONT> * @throws IllegalArgumentException if any parameters are invalid<a name="line.443"></a>
447 <FONT color="green">444</FONT> */<a name="line.444"></a>
448 <FONT color="green">445</FONT> public static double[] sample(UnivariateRealFunction f,<a name="line.445"></a>
449 <FONT color="green">446</FONT> double min, double max, int n)<a name="line.446"></a>
450 <FONT color="green">447</FONT> throws FunctionEvaluationException, IllegalArgumentException {<a name="line.447"></a>
451 <FONT color="green">448</FONT> <a name="line.448"></a>
452 <FONT color="green">449</FONT> if (n &lt;= 0) {<a name="line.449"></a>
453 <FONT color="green">450</FONT> throw MathRuntimeException.createIllegalArgumentException(<a name="line.450"></a>
454 <FONT color="green">451</FONT> "number of sample is not positive: {0}",<a name="line.451"></a>
455 <FONT color="green">452</FONT> n);<a name="line.452"></a>
456 <FONT color="green">453</FONT> }<a name="line.453"></a>
457 <FONT color="green">454</FONT> verifyInterval(min, max);<a name="line.454"></a>
458 <FONT color="green">455</FONT> <a name="line.455"></a>
459 <FONT color="green">456</FONT> double s[] = new double[n];<a name="line.456"></a>
460 <FONT color="green">457</FONT> double h = (max - min) / n;<a name="line.457"></a>
461 <FONT color="green">458</FONT> for (int i = 0; i &lt; n; i++) {<a name="line.458"></a>
462 <FONT color="green">459</FONT> s[i] = f.value(min + i * h);<a name="line.459"></a>
463 <FONT color="green">460</FONT> }<a name="line.460"></a>
464 <FONT color="green">461</FONT> return s;<a name="line.461"></a>
465 <FONT color="green">462</FONT> }<a name="line.462"></a>
466 <FONT color="green">463</FONT> <a name="line.463"></a>
467 <FONT color="green">464</FONT> /**<a name="line.464"></a>
468 <FONT color="green">465</FONT> * Multiply every component in the given real array by the<a name="line.465"></a>
469 <FONT color="green">466</FONT> * given real number. The change is made in place.<a name="line.466"></a>
470 <FONT color="green">467</FONT> *<a name="line.467"></a>
471 <FONT color="green">468</FONT> * @param f the real array to be scaled<a name="line.468"></a>
472 <FONT color="green">469</FONT> * @param d the real scaling coefficient<a name="line.469"></a>
473 <FONT color="green">470</FONT> * @return a reference to the scaled array<a name="line.470"></a>
474 <FONT color="green">471</FONT> */<a name="line.471"></a>
475 <FONT color="green">472</FONT> public static double[] scaleArray(double f[], double d) {<a name="line.472"></a>
476 <FONT color="green">473</FONT> for (int i = 0; i &lt; f.length; i++) {<a name="line.473"></a>
477 <FONT color="green">474</FONT> f[i] *= d;<a name="line.474"></a>
478 <FONT color="green">475</FONT> }<a name="line.475"></a>
479 <FONT color="green">476</FONT> return f;<a name="line.476"></a>
480 <FONT color="green">477</FONT> }<a name="line.477"></a>
481 <FONT color="green">478</FONT> <a name="line.478"></a>
482 <FONT color="green">479</FONT> /**<a name="line.479"></a>
483 <FONT color="green">480</FONT> * Multiply every component in the given complex array by the<a name="line.480"></a>
484 <FONT color="green">481</FONT> * given real number. The change is made in place.<a name="line.481"></a>
485 <FONT color="green">482</FONT> *<a name="line.482"></a>
486 <FONT color="green">483</FONT> * @param f the complex array to be scaled<a name="line.483"></a>
487 <FONT color="green">484</FONT> * @param d the real scaling coefficient<a name="line.484"></a>
488 <FONT color="green">485</FONT> * @return a reference to the scaled array<a name="line.485"></a>
489 <FONT color="green">486</FONT> */<a name="line.486"></a>
490 <FONT color="green">487</FONT> public static Complex[] scaleArray(Complex f[], double d) {<a name="line.487"></a>
491 <FONT color="green">488</FONT> for (int i = 0; i &lt; f.length; i++) {<a name="line.488"></a>
492 <FONT color="green">489</FONT> f[i] = new Complex(d * f[i].getReal(), d * f[i].getImaginary());<a name="line.489"></a>
493 <FONT color="green">490</FONT> }<a name="line.490"></a>
494 <FONT color="green">491</FONT> return f;<a name="line.491"></a>
495 <FONT color="green">492</FONT> }<a name="line.492"></a>
496 <FONT color="green">493</FONT> <a name="line.493"></a>
497 <FONT color="green">494</FONT> /**<a name="line.494"></a>
498 <FONT color="green">495</FONT> * Returns true if the argument is power of 2.<a name="line.495"></a>
499 <FONT color="green">496</FONT> *<a name="line.496"></a>
500 <FONT color="green">497</FONT> * @param n the number to test<a name="line.497"></a>
501 <FONT color="green">498</FONT> * @return true if the argument is power of 2<a name="line.498"></a>
502 <FONT color="green">499</FONT> */<a name="line.499"></a>
503 <FONT color="green">500</FONT> public static boolean isPowerOf2(long n) {<a name="line.500"></a>
504 <FONT color="green">501</FONT> return (n &gt; 0) &amp;&amp; ((n &amp; (n - 1)) == 0);<a name="line.501"></a>
505 <FONT color="green">502</FONT> }<a name="line.502"></a>
506 <FONT color="green">503</FONT> <a name="line.503"></a>
507 <FONT color="green">504</FONT> /**<a name="line.504"></a>
508 <FONT color="green">505</FONT> * Verifies that the data set has length of power of 2.<a name="line.505"></a>
509 <FONT color="green">506</FONT> *<a name="line.506"></a>
510 <FONT color="green">507</FONT> * @param d the data array<a name="line.507"></a>
511 <FONT color="green">508</FONT> * @throws IllegalArgumentException if array length is not power of 2<a name="line.508"></a>
512 <FONT color="green">509</FONT> */<a name="line.509"></a>
513 <FONT color="green">510</FONT> public static void verifyDataSet(double d[]) throws IllegalArgumentException {<a name="line.510"></a>
514 <FONT color="green">511</FONT> if (!isPowerOf2(d.length)) {<a name="line.511"></a>
515 <FONT color="green">512</FONT> throw MathRuntimeException.createIllegalArgumentException(<a name="line.512"></a>
516 <FONT color="green">513</FONT> NOT_POWER_OF_TWO_MESSAGE, d.length);<a name="line.513"></a>
517 <FONT color="green">514</FONT> }<a name="line.514"></a>
518 <FONT color="green">515</FONT> }<a name="line.515"></a>
519 <FONT color="green">516</FONT> <a name="line.516"></a>
520 <FONT color="green">517</FONT> /**<a name="line.517"></a>
521 <FONT color="green">518</FONT> * Verifies that the data set has length of power of 2.<a name="line.518"></a>
522 <FONT color="green">519</FONT> *<a name="line.519"></a>
523 <FONT color="green">520</FONT> * @param o the data array<a name="line.520"></a>
524 <FONT color="green">521</FONT> * @throws IllegalArgumentException if array length is not power of 2<a name="line.521"></a>
525 <FONT color="green">522</FONT> */<a name="line.522"></a>
526 <FONT color="green">523</FONT> public static void verifyDataSet(Object o[]) throws IllegalArgumentException {<a name="line.523"></a>
527 <FONT color="green">524</FONT> if (!isPowerOf2(o.length)) {<a name="line.524"></a>
528 <FONT color="green">525</FONT> throw MathRuntimeException.createIllegalArgumentException(<a name="line.525"></a>
529 <FONT color="green">526</FONT> NOT_POWER_OF_TWO_MESSAGE, o.length);<a name="line.526"></a>
530 <FONT color="green">527</FONT> }<a name="line.527"></a>
531 <FONT color="green">528</FONT> }<a name="line.528"></a>
532 <FONT color="green">529</FONT> <a name="line.529"></a>
533 <FONT color="green">530</FONT> /**<a name="line.530"></a>
534 <FONT color="green">531</FONT> * Verifies that the endpoints specify an interval.<a name="line.531"></a>
535 <FONT color="green">532</FONT> *<a name="line.532"></a>
536 <FONT color="green">533</FONT> * @param lower lower endpoint<a name="line.533"></a>
537 <FONT color="green">534</FONT> * @param upper upper endpoint<a name="line.534"></a>
538 <FONT color="green">535</FONT> * @throws IllegalArgumentException if not interval<a name="line.535"></a>
539 <FONT color="green">536</FONT> */<a name="line.536"></a>
540 <FONT color="green">537</FONT> public static void verifyInterval(double lower, double upper)<a name="line.537"></a>
541 <FONT color="green">538</FONT> throws IllegalArgumentException {<a name="line.538"></a>
542 <FONT color="green">539</FONT> <a name="line.539"></a>
543 <FONT color="green">540</FONT> if (lower &gt;= upper) {<a name="line.540"></a>
544 <FONT color="green">541</FONT> throw MathRuntimeException.createIllegalArgumentException(<a name="line.541"></a>
545 <FONT color="green">542</FONT> "endpoints do not specify an interval: [{0}, {1}]",<a name="line.542"></a>
546 <FONT color="green">543</FONT> lower, upper);<a name="line.543"></a>
547 <FONT color="green">544</FONT> }<a name="line.544"></a>
548 <FONT color="green">545</FONT> }<a name="line.545"></a>
549 <FONT color="green">546</FONT> <a name="line.546"></a>
550 <FONT color="green">547</FONT> /**<a name="line.547"></a>
551 <FONT color="green">548</FONT> * Performs a multi-dimensional Fourier transform on a given array.<a name="line.548"></a>
552 <FONT color="green">549</FONT> * Use {@link #inversetransform2(Complex[])} and<a name="line.549"></a>
553 <FONT color="green">550</FONT> * {@link #transform2(Complex[])} in a row-column implementation<a name="line.550"></a>
554 <FONT color="green">551</FONT> * in any number of dimensions with O(N&amp;times;log(N)) complexity with<a name="line.551"></a>
555 <FONT color="green">552</FONT> * N=n&lt;sub&gt;1&lt;/sub&gt;&amp;times;n&lt;sub&gt;2&lt;/sub&gt;&amp;times;n&lt;sub&gt;3&lt;/sub&gt;&amp;times;...&amp;times;n&lt;sub&gt;d&lt;/sub&gt;,<a name="line.552"></a>
556 <FONT color="green">553</FONT> * n&lt;sub&gt;x&lt;/sub&gt;=number of elements in dimension x,<a name="line.553"></a>
557 <FONT color="green">554</FONT> * and d=total number of dimensions.<a name="line.554"></a>
558 <FONT color="green">555</FONT> *<a name="line.555"></a>
559 <FONT color="green">556</FONT> * @param mdca Multi-Dimensional Complex Array id est Complex[][][][]<a name="line.556"></a>
560 <FONT color="green">557</FONT> * @param forward inverseTransform2 is preformed if this is false<a name="line.557"></a>
561 <FONT color="green">558</FONT> * @return transform of mdca as a Multi-Dimensional Complex Array id est Complex[][][][]<a name="line.558"></a>
562 <FONT color="green">559</FONT> * @throws IllegalArgumentException if any dimension is not a power of two<a name="line.559"></a>
563 <FONT color="green">560</FONT> */<a name="line.560"></a>
564 <FONT color="green">561</FONT> public Object mdfft(Object mdca, boolean forward)<a name="line.561"></a>
565 <FONT color="green">562</FONT> throws IllegalArgumentException {<a name="line.562"></a>
566 <FONT color="green">563</FONT> MultiDimensionalComplexMatrix mdcm = (MultiDimensionalComplexMatrix)<a name="line.563"></a>
567 <FONT color="green">564</FONT> new MultiDimensionalComplexMatrix(mdca).clone();<a name="line.564"></a>
568 <FONT color="green">565</FONT> int[] dimensionSize = mdcm.getDimensionSizes();<a name="line.565"></a>
569 <FONT color="green">566</FONT> //cycle through each dimension<a name="line.566"></a>
570 <FONT color="green">567</FONT> for (int i = 0; i &lt; dimensionSize.length; i++) {<a name="line.567"></a>
571 <FONT color="green">568</FONT> mdfft(mdcm, forward, i, new int[0]);<a name="line.568"></a>
572 <FONT color="green">569</FONT> }<a name="line.569"></a>
573 <FONT color="green">570</FONT> return mdcm.getArray();<a name="line.570"></a>
574 <FONT color="green">571</FONT> }<a name="line.571"></a>
575 <FONT color="green">572</FONT> <a name="line.572"></a>
576 <FONT color="green">573</FONT> /**<a name="line.573"></a>
577 <FONT color="green">574</FONT> * Performs one dimension of a multi-dimensional Fourier transform.<a name="line.574"></a>
578 <FONT color="green">575</FONT> *<a name="line.575"></a>
579 <FONT color="green">576</FONT> * @param mdcm input matrix<a name="line.576"></a>
580 <FONT color="green">577</FONT> * @param forward inverseTransform2 is preformed if this is false<a name="line.577"></a>
581 <FONT color="green">578</FONT> * @param d index of the dimension to process<a name="line.578"></a>
582 <FONT color="green">579</FONT> * @param subVector recursion subvector<a name="line.579"></a>
583 <FONT color="green">580</FONT> * @throws IllegalArgumentException if any dimension is not a power of two<a name="line.580"></a>
584 <FONT color="green">581</FONT> */<a name="line.581"></a>
585 <FONT color="green">582</FONT> private void mdfft(MultiDimensionalComplexMatrix mdcm, boolean forward,<a name="line.582"></a>
586 <FONT color="green">583</FONT> int d, int[] subVector)<a name="line.583"></a>
587 <FONT color="green">584</FONT> throws IllegalArgumentException {<a name="line.584"></a>
588 <FONT color="green">585</FONT> int[] dimensionSize = mdcm.getDimensionSizes();<a name="line.585"></a>
589 <FONT color="green">586</FONT> //if done<a name="line.586"></a>
590 <FONT color="green">587</FONT> if (subVector.length == dimensionSize.length) {<a name="line.587"></a>
591 <FONT color="green">588</FONT> Complex[] temp = new Complex[dimensionSize[d]];<a name="line.588"></a>
592 <FONT color="green">589</FONT> for (int i = 0; i &lt; dimensionSize[d]; i++) {<a name="line.589"></a>
593 <FONT color="green">590</FONT> //fft along dimension d<a name="line.590"></a>
594 <FONT color="green">591</FONT> subVector[d] = i;<a name="line.591"></a>
595 <FONT color="green">592</FONT> temp[i] = mdcm.get(subVector);<a name="line.592"></a>
596 <FONT color="green">593</FONT> }<a name="line.593"></a>
597 <FONT color="green">594</FONT> <a name="line.594"></a>
598 <FONT color="green">595</FONT> if (forward)<a name="line.595"></a>
599 <FONT color="green">596</FONT> temp = transform2(temp);<a name="line.596"></a>
600 <FONT color="green">597</FONT> else<a name="line.597"></a>
601 <FONT color="green">598</FONT> temp = inversetransform2(temp);<a name="line.598"></a>
602 <FONT color="green">599</FONT> <a name="line.599"></a>
603 <FONT color="green">600</FONT> for (int i = 0; i &lt; dimensionSize[d]; i++) {<a name="line.600"></a>
604 <FONT color="green">601</FONT> subVector[d] = i;<a name="line.601"></a>
605 <FONT color="green">602</FONT> mdcm.set(temp[i], subVector);<a name="line.602"></a>
606 <FONT color="green">603</FONT> }<a name="line.603"></a>
607 <FONT color="green">604</FONT> } else {<a name="line.604"></a>
608 <FONT color="green">605</FONT> int[] vector = new int[subVector.length + 1];<a name="line.605"></a>
609 <FONT color="green">606</FONT> System.arraycopy(subVector, 0, vector, 0, subVector.length);<a name="line.606"></a>
610 <FONT color="green">607</FONT> if (subVector.length == d) {<a name="line.607"></a>
611 <FONT color="green">608</FONT> //value is not important once the recursion is done.<a name="line.608"></a>
612 <FONT color="green">609</FONT> //then an fft will be applied along the dimension d.<a name="line.609"></a>
613 <FONT color="green">610</FONT> vector[d] = 0;<a name="line.610"></a>
614 <FONT color="green">611</FONT> mdfft(mdcm, forward, d, vector);<a name="line.611"></a>
615 <FONT color="green">612</FONT> } else {<a name="line.612"></a>
616 <FONT color="green">613</FONT> for (int i = 0; i &lt; dimensionSize[subVector.length]; i++) {<a name="line.613"></a>
617 <FONT color="green">614</FONT> vector[subVector.length] = i;<a name="line.614"></a>
618 <FONT color="green">615</FONT> //further split along the next dimension<a name="line.615"></a>
619 <FONT color="green">616</FONT> mdfft(mdcm, forward, d, vector);<a name="line.616"></a>
620 <FONT color="green">617</FONT> }<a name="line.617"></a>
621 <FONT color="green">618</FONT> }<a name="line.618"></a>
622 <FONT color="green">619</FONT> }<a name="line.619"></a>
623 <FONT color="green">620</FONT> return;<a name="line.620"></a>
624 <FONT color="green">621</FONT> }<a name="line.621"></a>
625 <FONT color="green">622</FONT> <a name="line.622"></a>
626 <FONT color="green">623</FONT> /**<a name="line.623"></a>
627 <FONT color="green">624</FONT> * Complex matrix implementation.<a name="line.624"></a>
628 <FONT color="green">625</FONT> * Not designed for synchronized access<a name="line.625"></a>
629 <FONT color="green">626</FONT> * may eventually be replaced by jsr-83 of the java community process<a name="line.626"></a>
630 <FONT color="green">627</FONT> * http://jcp.org/en/jsr/detail?id=83<a name="line.627"></a>
631 <FONT color="green">628</FONT> * may require additional exception throws for other basic requirements.<a name="line.628"></a>
632 <FONT color="green">629</FONT> */<a name="line.629"></a>
633 <FONT color="green">630</FONT> private static class MultiDimensionalComplexMatrix<a name="line.630"></a>
634 <FONT color="green">631</FONT> implements Cloneable {<a name="line.631"></a>
635 <FONT color="green">632</FONT> <a name="line.632"></a>
636 <FONT color="green">633</FONT> /** Size in all dimensions. */<a name="line.633"></a>
637 <FONT color="green">634</FONT> protected int[] dimensionSize;<a name="line.634"></a>
638 <FONT color="green">635</FONT> <a name="line.635"></a>
639 <FONT color="green">636</FONT> /** Storage array. */<a name="line.636"></a>
640 <FONT color="green">637</FONT> protected Object multiDimensionalComplexArray;<a name="line.637"></a>
641 <FONT color="green">638</FONT> <a name="line.638"></a>
642 <FONT color="green">639</FONT> /** Simple constructor.<a name="line.639"></a>
643 <FONT color="green">640</FONT> * @param multiDimensionalComplexArray array containing the matrix elements<a name="line.640"></a>
644 <FONT color="green">641</FONT> */<a name="line.641"></a>
645 <FONT color="green">642</FONT> public MultiDimensionalComplexMatrix(Object multiDimensionalComplexArray) {<a name="line.642"></a>
646 <FONT color="green">643</FONT> <a name="line.643"></a>
647 <FONT color="green">644</FONT> this.multiDimensionalComplexArray = multiDimensionalComplexArray;<a name="line.644"></a>
648 <FONT color="green">645</FONT> <a name="line.645"></a>
649 <FONT color="green">646</FONT> // count dimensions<a name="line.646"></a>
650 <FONT color="green">647</FONT> int numOfDimensions = 0;<a name="line.647"></a>
651 <FONT color="green">648</FONT> for (Object lastDimension = multiDimensionalComplexArray;<a name="line.648"></a>
652 <FONT color="green">649</FONT> lastDimension instanceof Object[];) {<a name="line.649"></a>
653 <FONT color="green">650</FONT> final Object[] array = (Object[]) lastDimension;<a name="line.650"></a>
654 <FONT color="green">651</FONT> numOfDimensions++;<a name="line.651"></a>
655 <FONT color="green">652</FONT> lastDimension = array[0];<a name="line.652"></a>
656 <FONT color="green">653</FONT> }<a name="line.653"></a>
657 <FONT color="green">654</FONT> <a name="line.654"></a>
658 <FONT color="green">655</FONT> // allocate array with exact count<a name="line.655"></a>
659 <FONT color="green">656</FONT> dimensionSize = new int[numOfDimensions];<a name="line.656"></a>
660 <FONT color="green">657</FONT> <a name="line.657"></a>
661 <FONT color="green">658</FONT> // fill array<a name="line.658"></a>
662 <FONT color="green">659</FONT> numOfDimensions = 0;<a name="line.659"></a>
663 <FONT color="green">660</FONT> for (Object lastDimension = multiDimensionalComplexArray;<a name="line.660"></a>
664 <FONT color="green">661</FONT> lastDimension instanceof Object[];) {<a name="line.661"></a>
665 <FONT color="green">662</FONT> final Object[] array = (Object[]) lastDimension;<a name="line.662"></a>
666 <FONT color="green">663</FONT> dimensionSize[numOfDimensions++] = array.length;<a name="line.663"></a>
667 <FONT color="green">664</FONT> lastDimension = array[0];<a name="line.664"></a>
668 <FONT color="green">665</FONT> }<a name="line.665"></a>
669 <FONT color="green">666</FONT> <a name="line.666"></a>
670 <FONT color="green">667</FONT> }<a name="line.667"></a>
671 <FONT color="green">668</FONT> <a name="line.668"></a>
672 <FONT color="green">669</FONT> /**<a name="line.669"></a>
673 <FONT color="green">670</FONT> * Get a matrix element.<a name="line.670"></a>
674 <FONT color="green">671</FONT> * @param vector indices of the element<a name="line.671"></a>
675 <FONT color="green">672</FONT> * @return matrix element<a name="line.672"></a>
676 <FONT color="green">673</FONT> * @exception IllegalArgumentException if dimensions do not match<a name="line.673"></a>
677 <FONT color="green">674</FONT> */<a name="line.674"></a>
678 <FONT color="green">675</FONT> public Complex get(int... vector)<a name="line.675"></a>
679 <FONT color="green">676</FONT> throws IllegalArgumentException {<a name="line.676"></a>
680 <FONT color="green">677</FONT> if (vector == null) {<a name="line.677"></a>
681 <FONT color="green">678</FONT> if (dimensionSize.length &gt; 0) {<a name="line.678"></a>
682 <FONT color="green">679</FONT> throw MathRuntimeException.createIllegalArgumentException(<a name="line.679"></a>
683 <FONT color="green">680</FONT> DIMENSION_MISMATCH_MESSAGE, 0, dimensionSize.length);<a name="line.680"></a>
684 <FONT color="green">681</FONT> }<a name="line.681"></a>
685 <FONT color="green">682</FONT> return null;<a name="line.682"></a>
686 <FONT color="green">683</FONT> }<a name="line.683"></a>
687 <FONT color="green">684</FONT> if (vector.length != dimensionSize.length) {<a name="line.684"></a>
688 <FONT color="green">685</FONT> throw MathRuntimeException.createIllegalArgumentException(<a name="line.685"></a>
689 <FONT color="green">686</FONT> DIMENSION_MISMATCH_MESSAGE, vector.length, dimensionSize.length);<a name="line.686"></a>
690 <FONT color="green">687</FONT> }<a name="line.687"></a>
691 <FONT color="green">688</FONT> <a name="line.688"></a>
692 <FONT color="green">689</FONT> Object lastDimension = multiDimensionalComplexArray;<a name="line.689"></a>
693 <FONT color="green">690</FONT> <a name="line.690"></a>
694 <FONT color="green">691</FONT> for (int i = 0; i &lt; dimensionSize.length; i++) {<a name="line.691"></a>
695 <FONT color="green">692</FONT> lastDimension = ((Object[]) lastDimension)[vector[i]];<a name="line.692"></a>
696 <FONT color="green">693</FONT> }<a name="line.693"></a>
697 <FONT color="green">694</FONT> return (Complex) lastDimension;<a name="line.694"></a>
698 <FONT color="green">695</FONT> }<a name="line.695"></a>
699 <FONT color="green">696</FONT> <a name="line.696"></a>
700 <FONT color="green">697</FONT> /**<a name="line.697"></a>
701 <FONT color="green">698</FONT> * Set a matrix element.<a name="line.698"></a>
702 <FONT color="green">699</FONT> * @param magnitude magnitude of the element<a name="line.699"></a>
703 <FONT color="green">700</FONT> * @param vector indices of the element<a name="line.700"></a>
704 <FONT color="green">701</FONT> * @return the previous value<a name="line.701"></a>
705 <FONT color="green">702</FONT> * @exception IllegalArgumentException if dimensions do not match<a name="line.702"></a>
706 <FONT color="green">703</FONT> */<a name="line.703"></a>
707 <FONT color="green">704</FONT> public Complex set(Complex magnitude, int... vector)<a name="line.704"></a>
708 <FONT color="green">705</FONT> throws IllegalArgumentException {<a name="line.705"></a>
709 <FONT color="green">706</FONT> if (vector == null) {<a name="line.706"></a>
710 <FONT color="green">707</FONT> if (dimensionSize.length &gt; 0) {<a name="line.707"></a>
711 <FONT color="green">708</FONT> throw MathRuntimeException.createIllegalArgumentException(<a name="line.708"></a>
712 <FONT color="green">709</FONT> DIMENSION_MISMATCH_MESSAGE, 0, dimensionSize.length);<a name="line.709"></a>
713 <FONT color="green">710</FONT> }<a name="line.710"></a>
714 <FONT color="green">711</FONT> return null;<a name="line.711"></a>
715 <FONT color="green">712</FONT> }<a name="line.712"></a>
716 <FONT color="green">713</FONT> if (vector.length != dimensionSize.length) {<a name="line.713"></a>
717 <FONT color="green">714</FONT> throw MathRuntimeException.createIllegalArgumentException(<a name="line.714"></a>
718 <FONT color="green">715</FONT> DIMENSION_MISMATCH_MESSAGE, vector.length,dimensionSize.length);<a name="line.715"></a>
719 <FONT color="green">716</FONT> }<a name="line.716"></a>
720 <FONT color="green">717</FONT> <a name="line.717"></a>
721 <FONT color="green">718</FONT> Object[] lastDimension = (Object[]) multiDimensionalComplexArray;<a name="line.718"></a>
722 <FONT color="green">719</FONT> for (int i = 0; i &lt; dimensionSize.length - 1; i++) {<a name="line.719"></a>
723 <FONT color="green">720</FONT> lastDimension = (Object[]) lastDimension[vector[i]];<a name="line.720"></a>
724 <FONT color="green">721</FONT> }<a name="line.721"></a>
725 <FONT color="green">722</FONT> <a name="line.722"></a>
726 <FONT color="green">723</FONT> Complex lastValue = (Complex) lastDimension[vector[dimensionSize.length - 1]];<a name="line.723"></a>
727 <FONT color="green">724</FONT> lastDimension[vector[dimensionSize.length - 1]] = magnitude;<a name="line.724"></a>
728 <FONT color="green">725</FONT> <a name="line.725"></a>
729 <FONT color="green">726</FONT> return lastValue;<a name="line.726"></a>
730 <FONT color="green">727</FONT> }<a name="line.727"></a>
731 <FONT color="green">728</FONT> <a name="line.728"></a>
732 <FONT color="green">729</FONT> /**<a name="line.729"></a>
733 <FONT color="green">730</FONT> * Get the size in all dimensions.<a name="line.730"></a>
734 <FONT color="green">731</FONT> * @return size in all dimensions<a name="line.731"></a>
735 <FONT color="green">732</FONT> */<a name="line.732"></a>
736 <FONT color="green">733</FONT> public int[] getDimensionSizes() {<a name="line.733"></a>
737 <FONT color="green">734</FONT> return dimensionSize.clone();<a name="line.734"></a>
738 <FONT color="green">735</FONT> }<a name="line.735"></a>
739 <FONT color="green">736</FONT> <a name="line.736"></a>
740 <FONT color="green">737</FONT> /**<a name="line.737"></a>
741 <FONT color="green">738</FONT> * Get the underlying storage array<a name="line.738"></a>
742 <FONT color="green">739</FONT> * @return underlying storage array<a name="line.739"></a>
743 <FONT color="green">740</FONT> */<a name="line.740"></a>
744 <FONT color="green">741</FONT> public Object getArray() {<a name="line.741"></a>
745 <FONT color="green">742</FONT> return multiDimensionalComplexArray;<a name="line.742"></a>
746 <FONT color="green">743</FONT> }<a name="line.743"></a>
747 <FONT color="green">744</FONT> <a name="line.744"></a>
748 <FONT color="green">745</FONT> /** {@inheritDoc} */<a name="line.745"></a>
749 <FONT color="green">746</FONT> @Override<a name="line.746"></a>
750 <FONT color="green">747</FONT> public Object clone() {<a name="line.747"></a>
751 <FONT color="green">748</FONT> MultiDimensionalComplexMatrix mdcm =<a name="line.748"></a>
752 <FONT color="green">749</FONT> new MultiDimensionalComplexMatrix(Array.newInstance(<a name="line.749"></a>
753 <FONT color="green">750</FONT> Complex.class, dimensionSize));<a name="line.750"></a>
754 <FONT color="green">751</FONT> clone(mdcm);<a name="line.751"></a>
755 <FONT color="green">752</FONT> return mdcm;<a name="line.752"></a>
756 <FONT color="green">753</FONT> }<a name="line.753"></a>
757 <FONT color="green">754</FONT> <a name="line.754"></a>
758 <FONT color="green">755</FONT> /**<a name="line.755"></a>
759 <FONT color="green">756</FONT> * Copy contents of current array into mdcm.<a name="line.756"></a>
760 <FONT color="green">757</FONT> * @param mdcm array where to copy data<a name="line.757"></a>
761 <FONT color="green">758</FONT> */<a name="line.758"></a>
762 <FONT color="green">759</FONT> private void clone(MultiDimensionalComplexMatrix mdcm) {<a name="line.759"></a>
763 <FONT color="green">760</FONT> int[] vector = new int[dimensionSize.length];<a name="line.760"></a>
764 <FONT color="green">761</FONT> int size = 1;<a name="line.761"></a>
765 <FONT color="green">762</FONT> for (int i = 0; i &lt; dimensionSize.length; i++) {<a name="line.762"></a>
766 <FONT color="green">763</FONT> size *= dimensionSize[i];<a name="line.763"></a>
767 <FONT color="green">764</FONT> }<a name="line.764"></a>
768 <FONT color="green">765</FONT> int[][] vectorList = new int[size][dimensionSize.length];<a name="line.765"></a>
769 <FONT color="green">766</FONT> for (int[] nextVector: vectorList) {<a name="line.766"></a>
770 <FONT color="green">767</FONT> System.arraycopy(vector, 0, nextVector, 0,<a name="line.767"></a>
771 <FONT color="green">768</FONT> dimensionSize.length);<a name="line.768"></a>
772 <FONT color="green">769</FONT> for (int i = 0; i &lt; dimensionSize.length; i++) {<a name="line.769"></a>
773 <FONT color="green">770</FONT> vector[i]++;<a name="line.770"></a>
774 <FONT color="green">771</FONT> if (vector[i] &lt; dimensionSize[i]) {<a name="line.771"></a>
775 <FONT color="green">772</FONT> break;<a name="line.772"></a>
776 <FONT color="green">773</FONT> } else {<a name="line.773"></a>
777 <FONT color="green">774</FONT> vector[i] = 0;<a name="line.774"></a>
778 <FONT color="green">775</FONT> }<a name="line.775"></a>
779 <FONT color="green">776</FONT> }<a name="line.776"></a>
780 <FONT color="green">777</FONT> }<a name="line.777"></a>
781 <FONT color="green">778</FONT> <a name="line.778"></a>
782 <FONT color="green">779</FONT> for (int[] nextVector: vectorList) {<a name="line.779"></a>
783 <FONT color="green">780</FONT> mdcm.set(get(nextVector), nextVector);<a name="line.780"></a>
784 <FONT color="green">781</FONT> }<a name="line.781"></a>
785 <FONT color="green">782</FONT> }<a name="line.782"></a>
786 <FONT color="green">783</FONT> }<a name="line.783"></a>
787 <FONT color="green">784</FONT> <a name="line.784"></a>
788 <FONT color="green">785</FONT> <a name="line.785"></a>
789 <FONT color="green">786</FONT> /** Computes the n&lt;sup&gt;th&lt;/sup&gt; roots of unity.<a name="line.786"></a>
790 <FONT color="green">787</FONT> * A cache of already computed values is maintained.<a name="line.787"></a>
791 <FONT color="green">788</FONT> */<a name="line.788"></a>
792 <FONT color="green">789</FONT> private static class RootsOfUnity implements Serializable {<a name="line.789"></a>
793 <FONT color="green">790</FONT> <a name="line.790"></a>
794 <FONT color="green">791</FONT> /** Serializable version id. */<a name="line.791"></a>
795 <FONT color="green">792</FONT> private static final long serialVersionUID = 6404784357747329667L;<a name="line.792"></a>
796 <FONT color="green">793</FONT> <a name="line.793"></a>
797 <FONT color="green">794</FONT> /** Number of roots of unity. */<a name="line.794"></a>
798 <FONT color="green">795</FONT> private int omegaCount;<a name="line.795"></a>
799 <FONT color="green">796</FONT> <a name="line.796"></a>
800 <FONT color="green">797</FONT> /** Real part of the roots. */<a name="line.797"></a>
801 <FONT color="green">798</FONT> private double[] omegaReal;<a name="line.798"></a>
802 <FONT color="green">799</FONT> <a name="line.799"></a>
803 <FONT color="green">800</FONT> /** Imaginary part of the roots for forward transform. */<a name="line.800"></a>
804 <FONT color="green">801</FONT> private double[] omegaImaginaryForward;<a name="line.801"></a>
805 <FONT color="green">802</FONT> <a name="line.802"></a>
806 <FONT color="green">803</FONT> /** Imaginary part of the roots for reverse transform. */<a name="line.803"></a>
807 <FONT color="green">804</FONT> private double[] omegaImaginaryInverse;<a name="line.804"></a>
808 <FONT color="green">805</FONT> <a name="line.805"></a>
809 <FONT color="green">806</FONT> /** Forward/reverse indicator. */<a name="line.806"></a>
810 <FONT color="green">807</FONT> private boolean isForward;<a name="line.807"></a>
811 <FONT color="green">808</FONT> <a name="line.808"></a>
812 <FONT color="green">809</FONT> /**<a name="line.809"></a>
813 <FONT color="green">810</FONT> * Build an engine for computing then &lt;sup&gt;th&lt;/sup&gt; roots of unity<a name="line.810"></a>
814 <FONT color="green">811</FONT> */<a name="line.811"></a>
815 <FONT color="green">812</FONT> public RootsOfUnity() {<a name="line.812"></a>
816 <FONT color="green">813</FONT> <a name="line.813"></a>
817 <FONT color="green">814</FONT> omegaCount = 0;<a name="line.814"></a>
818 <FONT color="green">815</FONT> omegaReal = null;<a name="line.815"></a>
819 <FONT color="green">816</FONT> omegaImaginaryForward = null;<a name="line.816"></a>
820 <FONT color="green">817</FONT> omegaImaginaryInverse = null;<a name="line.817"></a>
821 <FONT color="green">818</FONT> isForward = true;<a name="line.818"></a>
822 <FONT color="green">819</FONT> <a name="line.819"></a>
823 <FONT color="green">820</FONT> }<a name="line.820"></a>
824 <FONT color="green">821</FONT> <a name="line.821"></a>
825 <FONT color="green">822</FONT> /**<a name="line.822"></a>
826 <FONT color="green">823</FONT> * Check if computation has been done for forward or reverse transform.<a name="line.823"></a>
827 <FONT color="green">824</FONT> * @return true if computation has been done for forward transform<a name="line.824"></a>
828 <FONT color="green">825</FONT> * @throws IllegalStateException if no roots of unity have been computed yet<a name="line.825"></a>
829 <FONT color="green">826</FONT> */<a name="line.826"></a>
830 <FONT color="green">827</FONT> public synchronized boolean isForward() throws IllegalStateException {<a name="line.827"></a>
831 <FONT color="green">828</FONT> <a name="line.828"></a>
832 <FONT color="green">829</FONT> if (omegaCount == 0) {<a name="line.829"></a>
833 <FONT color="green">830</FONT> throw MathRuntimeException.createIllegalStateException(<a name="line.830"></a>
834 <FONT color="green">831</FONT> MISSING_ROOTS_OF_UNITY_MESSAGE);<a name="line.831"></a>
835 <FONT color="green">832</FONT> }<a name="line.832"></a>
836 <FONT color="green">833</FONT> return isForward;<a name="line.833"></a>
837 <FONT color="green">834</FONT> <a name="line.834"></a>
838 <FONT color="green">835</FONT> }<a name="line.835"></a>
839 <FONT color="green">836</FONT> <a name="line.836"></a>
840 <FONT color="green">837</FONT> /** Computes the n&lt;sup&gt;th&lt;/sup&gt; roots of unity.<a name="line.837"></a>
841 <FONT color="green">838</FONT> * &lt;p&gt;The computed omega[] = { 1, w, w&lt;sup&gt;2&lt;/sup&gt;, ... w&lt;sup&gt;(n-1)&lt;/sup&gt; } where<a name="line.838"></a>
842 <FONT color="green">839</FONT> * w = exp(-2 &amp;pi; i / n), i = &amp;sqrt;(-1).&lt;/p&gt;<a name="line.839"></a>
843 <FONT color="green">840</FONT> * &lt;p&gt;Note that n is positive for<a name="line.840"></a>
844 <FONT color="green">841</FONT> * forward transform and negative for inverse transform.&lt;/p&gt;<a name="line.841"></a>
845 <FONT color="green">842</FONT> * @param n number of roots of unity to compute,<a name="line.842"></a>
846 <FONT color="green">843</FONT> * positive for forward transform, negative for inverse transform<a name="line.843"></a>
847 <FONT color="green">844</FONT> * @throws IllegalArgumentException if n = 0<a name="line.844"></a>
848 <FONT color="green">845</FONT> */<a name="line.845"></a>
849 <FONT color="green">846</FONT> public synchronized void computeOmega(int n) throws IllegalArgumentException {<a name="line.846"></a>
850 <FONT color="green">847</FONT> <a name="line.847"></a>
851 <FONT color="green">848</FONT> if (n == 0) {<a name="line.848"></a>
852 <FONT color="green">849</FONT> throw MathRuntimeException.createIllegalArgumentException(<a name="line.849"></a>
853 <FONT color="green">850</FONT> "cannot compute 0-th root of unity, indefinite result");<a name="line.850"></a>
854 <FONT color="green">851</FONT> }<a name="line.851"></a>
855 <FONT color="green">852</FONT> <a name="line.852"></a>
856 <FONT color="green">853</FONT> isForward = n &gt; 0;<a name="line.853"></a>
857 <FONT color="green">854</FONT> <a name="line.854"></a>
858 <FONT color="green">855</FONT> // avoid repetitive calculations<a name="line.855"></a>
859 <FONT color="green">856</FONT> final int absN = Math.abs(n);<a name="line.856"></a>
860 <FONT color="green">857</FONT> <a name="line.857"></a>
861 <FONT color="green">858</FONT> if (absN == omegaCount) {<a name="line.858"></a>
862 <FONT color="green">859</FONT> return;<a name="line.859"></a>
863 <FONT color="green">860</FONT> }<a name="line.860"></a>
864 <FONT color="green">861</FONT> <a name="line.861"></a>
865 <FONT color="green">862</FONT> // calculate everything from scratch, for both forward and inverse versions<a name="line.862"></a>
866 <FONT color="green">863</FONT> final double t = 2.0 * Math.PI / absN;<a name="line.863"></a>
867 <FONT color="green">864</FONT> final double cosT = Math.cos(t);<a name="line.864"></a>
868 <FONT color="green">865</FONT> final double sinT = Math.sin(t);<a name="line.865"></a>
869 <FONT color="green">866</FONT> omegaReal = new double[absN];<a name="line.866"></a>
870 <FONT color="green">867</FONT> omegaImaginaryForward = new double[absN];<a name="line.867"></a>
871 <FONT color="green">868</FONT> omegaImaginaryInverse = new double[absN];<a name="line.868"></a>
872 <FONT color="green">869</FONT> omegaReal[0] = 1.0;<a name="line.869"></a>
873 <FONT color="green">870</FONT> omegaImaginaryForward[0] = 0.0;<a name="line.870"></a>
874 <FONT color="green">871</FONT> omegaImaginaryInverse[0] = 0.0;<a name="line.871"></a>
875 <FONT color="green">872</FONT> for (int i = 1; i &lt; absN; i++) {<a name="line.872"></a>
876 <FONT color="green">873</FONT> omegaReal[i] =<a name="line.873"></a>
877 <FONT color="green">874</FONT> omegaReal[i-1] * cosT + omegaImaginaryForward[i-1] * sinT;<a name="line.874"></a>
878 <FONT color="green">875</FONT> omegaImaginaryForward[i] =<a name="line.875"></a>
879 <FONT color="green">876</FONT> omegaImaginaryForward[i-1] * cosT - omegaReal[i-1] * sinT;<a name="line.876"></a>
880 <FONT color="green">877</FONT> omegaImaginaryInverse[i] = -omegaImaginaryForward[i];<a name="line.877"></a>
881 <FONT color="green">878</FONT> }<a name="line.878"></a>
882 <FONT color="green">879</FONT> omegaCount = absN;<a name="line.879"></a>
883 <FONT color="green">880</FONT> <a name="line.880"></a>
884 <FONT color="green">881</FONT> }<a name="line.881"></a>
885 <FONT color="green">882</FONT> <a name="line.882"></a>
886 <FONT color="green">883</FONT> /**<a name="line.883"></a>
887 <FONT color="green">884</FONT> * Get the real part of the k&lt;sup&gt;th&lt;/sup&gt; n&lt;sup&gt;th&lt;/sup&gt; root of unity<a name="line.884"></a>
888 <FONT color="green">885</FONT> * @param k index of the n&lt;sup&gt;th&lt;/sup&gt; root of unity<a name="line.885"></a>
889 <FONT color="green">886</FONT> * @return real part of the k&lt;sup&gt;th&lt;/sup&gt; n&lt;sup&gt;th&lt;/sup&gt; root of unity<a name="line.886"></a>
890 <FONT color="green">887</FONT> * @throws IllegalStateException if no roots of unity have been computed yet<a name="line.887"></a>
891 <FONT color="green">888</FONT> * @throws IllegalArgumentException if k is out of range<a name="line.888"></a>
892 <FONT color="green">889</FONT> */<a name="line.889"></a>
893 <FONT color="green">890</FONT> public synchronized double getOmegaReal(int k)<a name="line.890"></a>
894 <FONT color="green">891</FONT> throws IllegalStateException, IllegalArgumentException {<a name="line.891"></a>
895 <FONT color="green">892</FONT> <a name="line.892"></a>
896 <FONT color="green">893</FONT> if (omegaCount == 0) {<a name="line.893"></a>
897 <FONT color="green">894</FONT> throw MathRuntimeException.createIllegalStateException(<a name="line.894"></a>
898 <FONT color="green">895</FONT> MISSING_ROOTS_OF_UNITY_MESSAGE);<a name="line.895"></a>
899 <FONT color="green">896</FONT> }<a name="line.896"></a>
900 <FONT color="green">897</FONT> if ((k &lt; 0) || (k &gt;= omegaCount)) {<a name="line.897"></a>
901 <FONT color="green">898</FONT> throw MathRuntimeException.createIllegalArgumentException(<a name="line.898"></a>
902 <FONT color="green">899</FONT> OUT_OF_RANGE_ROOT_INDEX_MESSAGE, k, 0, omegaCount - 1);<a name="line.899"></a>
903 <FONT color="green">900</FONT> }<a name="line.900"></a>
904 <FONT color="green">901</FONT> <a name="line.901"></a>
905 <FONT color="green">902</FONT> return omegaReal[k];<a name="line.902"></a>
906 <FONT color="green">903</FONT> <a name="line.903"></a>
907 <FONT color="green">904</FONT> }<a name="line.904"></a>
908 <FONT color="green">905</FONT> <a name="line.905"></a>
909 <FONT color="green">906</FONT> /**<a name="line.906"></a>
910 <FONT color="green">907</FONT> * Get the imaginary part of the k&lt;sup&gt;th&lt;/sup&gt; n&lt;sup&gt;th&lt;/sup&gt; root of unity<a name="line.907"></a>
911 <FONT color="green">908</FONT> * @param k index of the n&lt;sup&gt;th&lt;/sup&gt; root of unity<a name="line.908"></a>
912 <FONT color="green">909</FONT> * @return imaginary part of the k&lt;sup&gt;th&lt;/sup&gt; n&lt;sup&gt;th&lt;/sup&gt; root of unity<a name="line.909"></a>
913 <FONT color="green">910</FONT> * @throws IllegalStateException if no roots of unity have been computed yet<a name="line.910"></a>
914 <FONT color="green">911</FONT> * @throws IllegalArgumentException if k is out of range<a name="line.911"></a>
915 <FONT color="green">912</FONT> */<a name="line.912"></a>
916 <FONT color="green">913</FONT> public synchronized double getOmegaImaginary(int k)<a name="line.913"></a>
917 <FONT color="green">914</FONT> throws IllegalStateException, IllegalArgumentException {<a name="line.914"></a>
918 <FONT color="green">915</FONT> <a name="line.915"></a>
919 <FONT color="green">916</FONT> if (omegaCount == 0) {<a name="line.916"></a>
920 <FONT color="green">917</FONT> throw MathRuntimeException.createIllegalStateException(<a name="line.917"></a>
921 <FONT color="green">918</FONT> MISSING_ROOTS_OF_UNITY_MESSAGE);<a name="line.918"></a>
922 <FONT color="green">919</FONT> }<a name="line.919"></a>
923 <FONT color="green">920</FONT> if ((k &lt; 0) || (k &gt;= omegaCount)) {<a name="line.920"></a>
924 <FONT color="green">921</FONT> throw MathRuntimeException.createIllegalArgumentException(<a name="line.921"></a>
925 <FONT color="green">922</FONT> OUT_OF_RANGE_ROOT_INDEX_MESSAGE, k, 0, omegaCount - 1);<a name="line.922"></a>
926 <FONT color="green">923</FONT> }<a name="line.923"></a>
927 <FONT color="green">924</FONT> <a name="line.924"></a>
928 <FONT color="green">925</FONT> return isForward ? omegaImaginaryForward[k] : omegaImaginaryInverse[k];<a name="line.925"></a>
929 <FONT color="green">926</FONT> <a name="line.926"></a>
930 <FONT color="green">927</FONT> }<a name="line.927"></a>
931 <FONT color="green">928</FONT> <a name="line.928"></a>
932 <FONT color="green">929</FONT> }<a name="line.929"></a>
933 <FONT color="green">930</FONT> <a name="line.930"></a>
934 <FONT color="green">931</FONT> }<a name="line.931"></a>
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995 </PRE>
996 </BODY>
997 </HTML>