1: package fourier;
2:
3: import java.util.Arrays;
4:
5: /**
6: * 離散フーリエ1次元変換(discrete fourier 1D transformation)。
7: * サンプル数(離散データ数)が2のN乗のときにはCooley-Tukeyのアルゴリズム(蝶(Butterfly)演算)を用いて高速化する。
8: * J.W.Cooley, J.W.Tukey: An Algorithm for the Machine Calculation of Complex Fourier Series, Mathematics of Computation, Vol.19 1965.
9: */
10: public class DiscreteFourier1dTransformation extends DiscreteFourierTransformation
11: {
12: /**
13: * 元データの実部を保持するフィールド。
14: */
15: protected double[] sourceRealPart;
16:
17: /**
18: * 元データの虚部を保持するフィールド。
19: */
20: protected double[] sourceImaginaryPart;
21:
22: /**
23: * フーリエ変換を施した実部を保持するフィールド。
24: */
25: protected double[] destinationRealPart;
26:
27: /**
28: * フーリエ変換を施した虚部を保持するフィールド。
29: */
30: protected double[] destinationImaginaryPart;
31:
32:
33: /**
34: * 元データの実部だけを指定して離散フーリエ1次元変換のインスタンスを作るコンストラクタ。
35: */
36: public DiscreteFourier1dTransformation(double[] sourceRealData)
37: {
38: super();
39: int size = sourceRealData.length;
40: double[] sourceImaginaryData = new double[size];
41: Arrays.fill(sourceImaginaryData, 0.0d);
42: this.initialize(sourceRealData, sourceImaginaryData);
43: }
44:
45: /**
46: * 元データの実部と虚部を指定して離散フーリエ1次元変換のインスタンスを作るコンストラクタ。
47: */
48: public DiscreteFourier1dTransformation(double[] sourceRealData, double[] sourceImaginaryData)
49: {
50: super();
51: this.initialize(sourceRealData, sourceImaginaryData);
52: }
53:
54: /**
55: * 離散フーリエ変換の結果を水に流す。
56: */
57: protected void flush()
58: {
59: destinationRealPart = null;
60: destinationImaginaryPart = null;
61: return;
62: }
63:
64: /**
65: * 離散フーリエ変換(順変換)を施した虚部を応答する。
66: */
67: public double[] imaginaryPart()
68: {
69: if (destinationImaginaryPart == null) { this.transform(); }
70: return destinationImaginaryPart;
71: }
72:
73: /**
74: * 離散フーリエ変換を初期化する。
75: */
76: protected void initialize()
77: {
78: sourceRealPart = null;
79: sourceImaginaryPart = null;
80: this.flush();
81: return;
82: }
83:
84: /**
85: * 離散フーリエ変換を初期化する。
86: */
87: protected void initialize(double[] sourceRealData, double[] sourceImaginaryData)
88: {
89: sourceRealPart = sourceRealData;
90: sourceImaginaryPart = sourceImaginaryData;
91: this.flush();
92: return;
93: }
94:
95: /**
96: * 逆離散フーリエ変換(逆変換)を施した虚部を応答する。
97: */
98: public double[] inverseImaginaryPart()
99: {
100: if (destinationImaginaryPart == null) { this.inverseTransform(); }
101: return destinationImaginaryPart;
102: }
103:
104: /**
105: * 逆離散フーリエ変換(逆変換)を施した実部を応答する。
106: */
107: public double[] inverseRealPart()
108: {
109: if (destinationRealPart == null) { this.inverseTransform(); }
110: return destinationRealPart;
111: }
112:
113: /**
114: * 逆離散フーリエ変換(逆変換)を施す。
115: */
116: public DiscreteFourier1dTransformation inverseTransform()
117: {
118: return this.transform(-1);
119: }
120:
121: /**
122: * ある整数(n)が2のN乗であるかどうかを応答する。
123: */
124: protected boolean isPowerOfTwo(int n)
125: {
126: return this.powerOfTwo(n) > 0;
127: }
128:
129: /**
130: * 配列(anArray)の対数を取って応答する。
131: */
132: public double[] logarithmicArray(double[] anArray)
133: {
134: int n = anArray.length;
135: double[] logarithmicArray = new double[n];
136: for (int i = 0; i < n; i++)
137: {
138: double aValue = anArray[i];
139: aValue = Math.log(aValue);
140: logarithmicArray[i] = aValue;
141: }
142: return logarithmicArray;
143: }
144:
145: /**
146: * 配列(anArray)を正規化して応答する。
147: */
148: public double[] normalizedArray(double[] anArray)
149: {
150: double maximum = 0.0d;
151: int n = anArray.length;
152: for (int i = 0; i < n; i++)
153: {
154: maximum = Math.max(anArray[i], maximum);
155: }
156: double[] normalizedArray = new double[n];
157: for (int i = 0; i < n; i++)
158: {
159: double aValue = anArray[i] / maximum;
160: aValue = Math.max(aValue, 0.0d);
161: aValue = Math.min(aValue, 1.0d);
162: normalizedArray[i] = aValue;
163: }
164: return normalizedArray;
165: }
166:
167: /**
168: * 離散フーリエ変換を施した結果の正規化対数パワースペクトルを応答する。
169: */
170: public double[] normalizedLogarithmicPowerSpectrum()
171: {
172: return normalizedArray(logarithmicArray(this.powerSpectrum()));
173: }
174:
175: /**
176: * 離散フーリエ変換を施した結果の正規化対数スペクトルを応答する。
177: */
178: public double[] normalizedLogarithmicSpectrum()
179: {
180: return normalizedArray(logarithmicArray(this.spectrum()));
181: }
182:
183: /**
184: * 離散フーリエ変換を施した結果の正規化パワースペクトルを応答する。
185: */
186: public double[] normalizedPowerSpectrum()
187: {
188: return normalizedArray(this.powerSpectrum());
189: }
190:
191: /**
192: * 離散フーリエ変換を施した結果の正規化スペクトルを応答する。
193: */
194: public double[] normalizedSpectrum()
195: {
196: return normalizedArray(this.spectrum());
197: }
198:
199: /**
200: * ある整数(n)が2のN乗であるとき、そのNを応答する。それ以外のときは、0を応答する。
201: */
202: protected int powerOfTwo(int n)
203: {
204: int powerOfTwo = 1;
205: for (int nth = 1; powerOfTwo < n; nth++)
206: {
207: powerOfTwo = 1 << nth;
208: if (powerOfTwo == n) { return nth; }
209: }
210: return 0;
211: }
212:
213: /**
214: * 離散フーリエ変換を施した結果のパワースペクトルを応答する。
215: */
216: public double[] powerSpectrum()
217: {
218: double[] realValues = this.realPart();
219: double[] imaginaryValues = this.imaginaryPart();
220: int n = realValues.length;
221: double[] powerSpectrum = new double[n];
222: for (int i = 0; i < n; i++)
223: {
224: double realValue = realValues[i];
225: double imaginaryValue = imaginaryValues[i];
226: double aValue = (realValue * realValue) + (imaginaryValue * imaginaryValue);
227: powerSpectrum[i] = aValue;
228: }
229: return powerSpectrum;
230: }
231:
232: /**
233: * 離散フーリエ変換(順変換)を施した実部を応答する。
234: */
235: public double[] realPart()
236: {
237: if (destinationRealPart == null) { this.transform(); }
238: return destinationRealPart;
239: }
240:
241: /**
242: * 離散フーリエ変換を施した結果のスペクトルを応答する。
243: */
244: public double[] spectrum()
245: {
246: return squareRootArray(this.powerSpectrum());
247: }
248:
249: /**
250: * 配列(anArray)の平方根を取って応答する。
251: */
252: public double[] squareRootArray(double[] anArray)
253: {
254: int n = anArray.length;
255: double[] squareRootArray = new double[n];
256: for (int i = 0; i < n; i++)
257: {
258: double aValue = anArray[i];
259: aValue = Math.sqrt(aValue);
260: squareRootArray[i] = aValue;
261: }
262: return squareRootArray;
263: }
264:
265: /**
266: * 離散フーリエ変換(順変換)を施す。
267: */
268: public DiscreteFourier1dTransformation transform()
269: {
270: return this.transform(1);
271: }
272:
273: /**
274: * 離散フーリエ変換を施す。flagが1の時は順変換。flagが-1の時は逆変換。
275: */
276: protected DiscreteFourier1dTransformation transform(int flag)
277: {
278: int n = sourceRealPart.length;
279: if (this.isPowerOfTwo(n)) { return this.turboTransform(flag); }
280: double[] x = this.sourceRealPart;
281: double[] y = this.sourceImaginaryPart;
282: double[] c = new double[n];
283: double[] s = new double[n];
284: double q = 2.0d * Math.PI / (double)n * (double)flag;
285: for (int i = 0; i < n; i++)
286: {
287: c[i] = 0.0d;
288: s[i] = 0.0d;
289: for (int j = 0; j < n; j++)
290: {
291: double xj = x[j];
292: double yj = y[j];
293: double qij = q * i * j;
294: double cos = Math.cos(qij);
295: double sin = Math.sin(qij);
296: c[i] += (xj * cos) - (yj * sin);
297: s[i] += (xj * sin) + (yj * cos);
298: }
299: if (flag == -1)
300: {
301: c[i] /= (double)n;
302: s[i] /= (double)n;
303: }
304: }
305: destinationRealPart = c;
306: destinationImaginaryPart = s;
307: return this;
308: }
309:
310: /**
311: * 高速離散フーリエ変換を施す。flagが1の時は順変換。flagが-1の時は逆変換。
312: * Cooley-Tukeyのアルゴリズム(蝶(Butterfly)演算)で高速化する。
313: * J.W.Cooley, J.W.Tukey: An Algorithm for the Machine Calculation of Complex Fourier Series, Mathematics of Computation, Vol.19 1965.
314: */
315: protected DiscreteFourier1dTransformation turboTransform(int flag)
316: {
317: int n = sourceRealPart.length;
318: int power = this.powerOfTwo(n);
319: double[] x = new double[n];
320: double[] y = new double[n];
321: System.arraycopy(this.sourceRealPart, 0, x, 0, n);
322: System.arraycopy(this.sourceImaginaryPart, 0, y, 0, n);
323: double[] c = new double[n];
324: double[] s = new double[n];
325: if (flag == -1)
326: {
327: for (int i = 0; i < n; i++)
328: {
329: x[i] /= (double)n;
330: y[i] /= (double)n;
331: }
332: }
333: this.turboTransformInitialize(x, y, power);
334: double unitAngle = 2.0d * Math.PI / (double)n;
335: int dft = 2;
336: for (int i = 0; i < power; i++)
337: {
338: int dftN = n / dft;
339: double stepAngle = unitAngle * (double)dftN;
340: for (int j = 0; j < dftN; j++)
341: {
342: double angle = 0.0d;
343: for (int k = 0; k < dft; k++)
344: {
345: int toIndex = j * dft + k;
346: int fromIndex1;
347: int fromIndex2;
348: int dftHalf = dft / 2;
349: if (k < dftHalf)
350: {
351: fromIndex1 = toIndex;
352: fromIndex2 = fromIndex1 + dftHalf;
353: }
354: else
355: {
356: fromIndex2 = toIndex;
357: fromIndex1 = toIndex - dftHalf;
358: }
359: double xr = x[fromIndex2];
360: double yi = y[fromIndex2];
361: double cos = Math.cos(angle);
362: double sin = Math.sin(angle);
363: c[toIndex] = x[fromIndex1] + (xr * cos) - ((double)flag * yi * sin);
364: s[toIndex] = y[fromIndex1] + ((double)flag * xr * sin) + (yi * cos);
365: angle += stepAngle;
366: }
367: }
368: for (int j = 0; j < n; j++)
369: {
370: x[j] = c[j];
371: y[j] = s[j];
372: }
373: dft *= 2;
374: }
375: destinationRealPart = x;
376: destinationImaginaryPart = y;
377: return this;
378: }
379:
380: /**
381: * 高速離散フーリエ変換のための置き換え。
382: * Cooley-Tukeyのアルゴリズムを適用するための準備。
383: */
384: private void turboTransformInitialize(double[] x, double[] y, int power)
385: {
386: int n = x.length;
387: double[] xr = new double[n];
388: double[] yi = new double[n];
389: int dftTimes = n;
390: for (int i = 0; i < power - 1; i++)
391: {
392: int toIndex = 0;
393: int offset = 0;
394: while (toIndex < n)
395: {
396: int fromIndex = 0;
397: while (fromIndex < dftTimes)
398: {
399: xr[toIndex] = x[fromIndex + offset];
400: yi[toIndex] = y[fromIndex + offset];
401: toIndex++;
402: fromIndex += 2;
403: if (fromIndex == dftTimes) { fromIndex = 1; }
404: }
405: offset += dftTimes;
406: }
407: for (int j = 0; j < n; j++)
408: {
409: x[j] = xr[j];
410: y[j] = yi[j];
411: }
412: dftTimes /= 2;
413: }
414: }
415: }
This document was generated by KSU.TextDoclet on 2012/05/28 at 06:08:28.