FFT_original.as

リビジョン 3165, 29.6 kB (コミッタ: keim, コミット時期: 3 年前)
SiON update to ver0.58, add FFT module, change voice architecture and many bug fixes

Line
1	//----------------------------------------------------------------------------------------------------
2	// Fast Fourier Transform module
3	// Ported and modified by keim.
4	// This soruce code is distributed under BSD-style license (see org.si.license.txt).
5	//
6	// Original code (written by c)
7	// The original source code is free licensed.
8	//----- < Following text is from original code's readme.txt. > -----
9	// Copyright(C) 1996-2001 Takuya OOURA
10	// email: ooura@mmm.t.u-tokyo.ac.jp
11	// download: http://momonga.t.u-tokyo.ac.jp/~ooura/fft.html
12	// You may use, copy, modify this code for any purpose and
13	// without fee. You may distribute this ORIGINAL package.
14	//----- < up to here > -----
15	// [NOTES] Now the download site is moved to http://www.kurims.kyoto-u.ac.jp/~ooura/fft.html.
16	// And the email address might be unavailable.
17	//----------------------------------------------------------------------------------------------------
18
19
20
21
22	package org.si.utils {
23	/** Fast Fourier Transform module */
24	public class FFT_original
25	{
26	// valiables
27	//------------------------------------------------------------
28	private var _length:int = 0;
29	private var _waveTable:Vector.<Number> = new Vector.<Number>();
30	private var _cosTable :Vector.<Number> = new Vector.<Number>();
31	private var _bitrvTemp:Vector.<int> = new Vector.<int>(256);
32
33
34
35
36	// constructor
37	//------------------------------------------------------------
38	/** constructor, specify calculating data length first. */
39	function FFT_original(length:int)
40	{
41	_initialize(length);
42	}
43
44
45
46
47	// calculation (original functions)
48	//------------------------------------------------------------
49	/** Complex Discrete Fourier Tranform
50	* @param isgn 1 for FFT, -1 for IFFT.
51	* @param src data to transform. The length must be same as you passed to constructor.
52	*/ public function cdft(isgn:int, src:Vector.<Number>) : void
53	{
54	if (isgn >= 0) {
55	_bitrv2(src, _length);
56	_cftfsub(src);
57	} else {
58	_bitrv2conj(src, _length);
59	_cftbsub(src);
60	}
61	}
62
63
64
65	/** Real Discrete Fourier Tranform
66	* @param isgn 1 for FFT, -1 for IFFT.
67	* @param src data to transform. The length must be same as you passed to constructor.
68	*/
69	public function rdft(isgn:int, src:Vector.<Number>) : void
70	{
71	var xi:Number;
72
73	if (isgn >= 0) {
74	_bitrv2(src, _length);
75	_cftfsub(src);
76	_rftfsub(src);
77	xi = src[0] - src[1];
78	src[0] += src[1];
79	src[1] = xi;
80	} else {
81	src[1] = 0.5 * (src[0] - src[1]);
82	src[0] -= src[1];
83	_rftbsub(src);
84	_bitrv2(src, _length);
85	_cftbsub(src);
86	}
87	}
88
89
90	/** Discrete Cosine Tranform
91	* @param isgn [ATTENTION] -1 for DCT, 1 for IDCT. Opposite from FFT.
92	* @param src data to transform. The length must be same as you passed to constructor.
93	*/
94	public function ddct(isgn:int, src:Vector.<Number>) : void
95	{
96	var j:int, xr:Number;
97
98	if (isgn < 0) {
99	xr = src[_length - 1];
100	for (j = _length - 2; j >= 2; j -= 2) {
101	src[j+1] = src[j] - src[j - 1];
102	src[j] += src[j - 1];
103	}
104	src[1] = src[0] - xr;
105	src[0] += xr;
106	_rftbsub(src);
107	_bitrv2(src, _length);
108	_cftbsub(src);
109	_dctsub(src);
110	} else {
111	_dctsub(src);
112	_bitrv2(src, _length);
113	_cftfsub(src);
114	_rftfsub(src);
115	xr = src[0] - src[1];
116	src[0] += src[1];
117	for (j = 2; j < _length; j += 2) {
118	src[j - 1] = src[j] - src[j+1];
119	src[j] += src[j+1];
120	}
121	src[_length - 1] = xr;
122	}
123	}
124
125
126	/** Discrete Sine Tranform
127	* @param isgn [ATTENTION] -1 for DST, 1 for IDST. Opposite from FFT.
128	* @param src data to transform. The length must be same as you passed to constructor.
129	*/
130	public function ddst(isgn:int, src:Vector.<Number>) : void
131	{
132	var j:int, xr:Number;
133
134	if (isgn < 0) {
135	xr = src[_length - 1];
136	for (j = _length - 2; j >= 2; j -= 2) {
137	src[j+1] = -src[j] - src[j - 1];
138	src[j] -= src[j - 1];
139	}
140	src[1] = src[0] + xr;
141	src[0] -= xr;
142	_rftbsub(src);
143	_bitrv2(src, _length);
144	_cftbsub(src);
145	_dstsub(src);
146	} else {
147	_dstsub(src);
148	_bitrv2(src, _length);
149	_cftfsub(src);
150	_rftfsub(src);
151	xr = src[0] - src[1];
152	src[0] += src[1];
153	for (j = 2; j < _length; j += 2) {
154	src[j - 1] = -src[j] - src[j+1];
155	src[j] -= src[j+1];
156	}
157	src[_length - 1] = -xr;
158	}
159	}
160
161
162
163
164	// internal function
165	//------------------------------------------------------------
166	// initializer
167	private function _initialize(len:int) : void
168	{
169	for (_length=8; _length<len;) _length<<=1;
170	_waveTable.length = _length >> 2;
171	_cosTable.length = _length;
172	var i:int, imax:int = _length >> 3, tlen:int = _waveTable.length,
173	dt:Number = 6.283185307179586 / _length;
174
175	_waveTable[0] = 1;
176	_waveTable[1] = 0;
177	_waveTable[imax+1] = _waveTable[imax] = Math.cos(0.7853981633974483);
178	for (i=2; i<imax; i+=2) {
179	_waveTable[tlen-i+1] = _waveTable[i] = Math.cos(i*dt);
180	_waveTable[tlen-i] = _waveTable[i+1] = Math.sin(i*dt);
181	}
182	_bitrv2(_waveTable, tlen);
183
184	imax = _cosTable.length;
185	dt = 1.5707963267948965 / imax;
186	for (i=0; i<imax; i++) _cosTable[i] = Math.cos(idt) 0.5;
187	}
188
189
190	// bit reverse
191	private function _bitrv2(src:Vector.<Number>, srclen:int) : void
192	{
193	var j:int, j1:int, k:int, k1:int,
194	xr:Number, xi:Number, yr:Number, yi:Number;
195
196	_bitrvTemp[0] = 0;
197	var l:int = srclen, m:int = 1;
198	while ((m << 3) < l) {
199	l >>= 1;
200	for (j = 0; j < m; j++) _bitrvTemp[m + j] = _bitrvTemp[j] + l;
201	m <<= 1;
202	}
203	var m2:int = m * 2;
204
205	if ((m << 3) == l) {
206	for (k = 0; k < m; k++) {
207	for (j = 0; j < k; j++) {
208	j1 = j + j + _bitrvTemp[k];
209	k1 = k + k + _bitrvTemp[j];
210	xr = src[j1];
211	xi = src[j1+1];
212	yr = src[k1];
213	yi = src[k1+1];
214	src[j1] = yr;
215	src[j1+1] = yi;
216	src[k1] = xr;
217	src[k1+1] = xi;
218	j1 += m2;
219	k1 += m2 + m2;
220	xr = src[j1];
221	xi = src[j1+1];
222	yr = src[k1];
223	yi = src[k1+1];
224	src[j1] = yr;
225	src[j1+1] = yi;
226	src[k1] = xr;
227	src[k1+1] = xi;
228	j1 += m2;
229	k1 -= m2;
230	xr = src[j1];
231	xi = src[j1+1];
232	yr = src[k1];
233	yi = src[k1+1];
234	src[j1] = yr;
235	src[j1+1] = yi;
236	src[k1] = xr;
237	src[k1+1] = xi;
238	j1 += m2;
239	k1 += m2 + m2;
240	xr = src[j1];
241	xi = src[j1+1];
242	yr = src[k1];
243	yi = src[k1+1];
244	src[j1] = yr;
245	src[j1+1] = yi;
246	src[k1] = xr;
247	src[k1+1] = xi;
248	}
249	j1 = k + k + m2 + _bitrvTemp[k];
250	k1 = j1 + m2;
251	xr = src[j1];
252	xi = src[j1+1];
253	yr = src[k1];
254	yi = src[k1+1];
255	src[j1] = yr;
256	src[j1+1] = yi;
257	src[k1] = xr;
258	src[k1+1] = xi;
259	}
260	} else {
261	for (k = 1; k < m; k++) {
262	for (j = 0; j < k; j++) {
263	j1 = j + j + _bitrvTemp[k];
264	k1 = k + k + _bitrvTemp[j];
265	xr = src[j1];
266	xi = src[j1+1];
267	yr = src[k1];
268	yi = src[k1+1];
269	src[j1] = yr;
270	src[j1+1] = yi;
271	src[k1] = xr;
272	src[k1+1] = xi;
273	j1 += m2;
274	k1 += m2;
275	xr = src[j1];
276	xi = src[j1+1];
277	yr = src[k1];
278	yi = src[k1+1];
279	src[j1] = yr;
280	src[j1+1] = yi;
281	src[k1] = xr;
282	src[k1+1] = xi;
283	}
284	}
285	}
286	}
287
288
289	// bit reverse (conjugation)
290	private function _bitrv2conj(src:Vector.<Number>, srclen:int) : void
291	{
292	var j:int, j1:int, k:int, k1:int,
293	xr:Number, xi:Number, yr:Number, yi:Number;
294
295	_bitrvTemp[0] = 0;
296	var l:int = srclen, m:int = 1;
297	while ((m << 3) < l) {
298	l >>= 1;
299	for (j = 0; j < m; j++) _bitrvTemp[m + j] = _bitrvTemp[j] + l;
300	m <<= 1;
301	}
302	var m2:int = m << 1;
303
304	if ((m << 3) == l) {
305	for (k = 0; k < m; k++) {
306	for (j = 0; j < k; j++) {
307	j1 = j + j + _bitrvTemp[k];
308	k1 = k + k + _bitrvTemp[j];
309	xr = src[j1];
310	xi = -src[j1+1];
311	yr = src[k1];
312	yi = -src[k1+1];
313	src[j1] = yr;
314	src[j1+1] = yi;
315	src[k1] = xr;
316	src[k1+1] = xi;
317	j1 += m2;
318	k1 += m2 + m2;
319	xr = src[j1];
320	xi = -src[j1+1];
321	yr = src[k1];
322	yi = -src[k1+1];
323	src[j1] = yr;
324	src[j1+1] = yi;
325	src[k1] = xr;
326	src[k1+1] = xi;
327	j1 += m2;
328	k1 -= m2;
329	xr = src[j1];
330	xi = -src[j1+1];
331	yr = src[k1];
332	yi = -src[k1+1];
333	src[j1] = yr;
334	src[j1+1] = yi;
335	src[k1] = xr;
336	src[k1+1] = xi;
337	j1 += m2;
338	k1 += m2 + m2;
339	xr = src[j1];
340	xi = -src[j1+1];
341	yr = src[k1];
342	yi = -src[k1+1];
343	src[j1] = yr;
344	src[j1+1] = yi;
345	src[k1] = xr;
346	src[k1+1] = xi;
347	}
348	k1 = k + k + _bitrvTemp[k];
349	src[k1+1] = -src[k1+1];
350	j1 = k1 + m2;
351	k1 = j1 + m2;
352	xr = src[j1];
353	xi = -src[j1+1];
354	yr = src[k1];
355	yi = -src[k1+1];
356	src[j1] = yr;
357	src[j1+1] = yi;
358	src[k1] = xr;
359	src[k1+1] = xi;
360	k1 += m2;
361	src[k1+1] = -src[k1+1];
362	}
363	} else {
364	src[1] = -src[1];
365	src[m2+1] = -src[m2+1];
366	for (k = 1; k < m; k++) {
367	for (j = 0; j < k; j++) {
368	j1 = j + j + _bitrvTemp[k];
369	k1 = k + k + _bitrvTemp[j];
370	xr = src[j1];
371	xi = -src[j1+1];
372	yr = src[k1];
373	yi = -src[k1+1];
374	src[j1] = yr;
375	src[j1+1] = yi;
376	src[k1] = xr;
377	src[k1+1] = xi;
378	j1 += m2;
379	k1 += m2;
380	xr = src[j1];
381	xi = -src[j1+1];
382	yr = src[k1];
383	yi = -src[k1+1];
384	src[j1] = yr;
385	src[j1+1] = yi;
386	src[k1] = xr;
387	src[k1+1] = xi;
388	}
389	k1 = k + k + _bitrvTemp[k];
390	src[k1+1] = -src[k1+1];
391	src[k1 + m2+1] = -src[k1 + m2+1];
392	}
393	}
394	}
395
396
397
398
399	// sub routines
400	//------------------------------------------------------------
401	private function _cftfsub(src:Vector.<Number>) : void
402	{
403	var j0:int, j1:int, j2:int, j3:int, l:int,
404	x0r:Number, x1r:Number, x2r:Number, x3r:Number,
405	x0i:Number, x1i:Number, x2i:Number, x3i:Number;
406
407	_cft1st(src);
408	l = 8;
409	while ((l << 2) < _length) {
410	_cftmdl(src, l);
411	l <<= 2;
412	}
413
414	if ((l << 2) == _length) {
415	for (j0 = 0; j0 < l; j0 += 2) {
416	j1 = j0 + l;
417	j2 = j1 + l;
418	j3 = j2 + l;
419	x0r = src[j0] + src[j1];
420	x0i = src[j0+1] + src[j1+1];
421	x1r = src[j0] - src[j1];
422	x1i = src[j0+1] - src[j1+1];
423	x2r = src[j2] + src[j3];
424	x2i = src[j2+1] + src[j3+1];
425	x3r = src[j2] - src[j3];
426	x3i = src[j2+1] - src[j3+1];
427	src[j0] = x0r + x2r;
428	src[j0+1] = x0i + x2i;
429	src[j2] = x0r - x2r;
430	src[j2+1] = x0i - x2i;
431	src[j1] = x1r - x3i;
432	src[j1+1] = x1i + x3r;
433	src[j3] = x1r + x3i;
434	src[j3+1] = x1i - x3r;
435	}
436	} else {
437	for (j0 = 0; j0 < l; j0 += 2) {
438	j1 = j0 + l;
439	x0r = src[j0] - src[j1];
440	x0i = src[j0+1] - src[j1+1];
441	src[j0] += src[j1];
442	src[j0+1] += src[j1+1];
443	src[j1] = x0r;
444	src[j1+1] = x0i;
445	}
446	}
447	}
448
449
450	private function _cftbsub(src:Vector.<Number>) : void
451	{
452	var j0:int, j1:int, j2:int, j3:int, l:int,
453	x0r:Number, x1r:Number, x2r:Number, x3r:Number,
454	x0i:Number, x1i:Number, x2i:Number, x3i:Number;
455
456	_cft1st(src);
457	l = 8;
458	while ((l << 2) < _length) {
459	_cftmdl(src, l);
460	l <<= 2;
461	}
462
463	if ((l << 2) == _length) {
464	for (j0 = 0; j0 < l; j0 += 2) {
465	j1 = j0 + l;
466	j2 = j1 + l;
467	j3 = j2 + l;
468	x0r = src[j0] + src[j1];
469	x0i = -src[j0+1] - src[j1+1];
470	x1r = src[j0] - src[j1];
471	x1i = -src[j0+1] + src[j1+1];
472	x2r = src[j2] + src[j3];
473	x2i = src[j2+1] + src[j3+1];
474	x3r = src[j2] - src[j3];
475	x3i = src[j2+1] - src[j3+1];
476	src[j0] = x0r + x2r;
477	src[j0+1] = x0i - x2i;
478	src[j2] = x0r - x2r;
479	src[j2+1] = x0i + x2i;
480	src[j1] = x1r - x3i;
481	src[j1+1] = x1i - x3r;
482	src[j3] = x1r + x3i;
483	src[j3+1] = x1i + x3r;
484	}
485	} else {
486	for (j0 = 0; j0 < l; j0 += 2) {
487	j1 = j0 + l;
488	x0r = src[j0] - src[j1];
489	x0i = -src[j0+1] + src[j1+1];
490	src[j0] += src[j1];
491	src[j0+1] = -src[j0+1] - src[j1+1];
492	src[j1] = x0r;
493	src[j1+1] = x0i;
494	}
495	}
496	}
497
498
499	private function _cft1st(src:Vector.<Number>) : void
500	{
501	var j:int, k1:int, k2:int,
502	wk1r:Number, wk2r:Number, wk3r:Number, x0r:Number, x1r:Number, x2r:Number, x3r:Number,
503	wk1i:Number, wk2i:Number, wk3i:Number, x0i:Number, x1i:Number, x2i:Number, x3i:Number;
504
505	x0r = src[0] + src[2];
506	x0i = src[1] + src[3];
507	x1r = src[0] - src[2];
508	x1i = src[1] - src[3];
509	x2r = src[4] + src[6];
510	x2i = src[5] + src[7];
511	x3r = src[4] - src[6];
512	x3i = src[5] - src[7];
513	src[0] = x0r + x2r;
514	src[1] = x0i + x2i;
515	src[4] = x0r - x2r;
516	src[5] = x0i - x2i;
517	src[2] = x1r - x3i;
518	src[3] = x1i + x3r;
519	src[6] = x1r + x3i;
520	src[7] = x1i - x3r;
521	wk1r = _waveTable[2];
522	x0r = src[8] + src[10];
523	x0i = src[9] + src[11];
524	x1r = src[8] - src[10];
525	x1i = src[9] - src[11];
526	x2r = src[12] + src[14];
527	x2i = src[13] + src[15];
528	x3r = src[12] - src[14];
529	x3i = src[13] - src[15];
530	src[8] = x0r + x2r;
531	src[9] = x0i + x2i;
532	src[12] = x2i - x0i;
533	src[13] = x0r - x2r;
534	x0r = x1r - x3i;
535	x0i = x1i + x3r;
536	src[10] = wk1r * (x0r - x0i);
537	src[11] = wk1r * (x0r + x0i);
538	x0r = x3i + x1r;
539	x0i = x3r - x1i;
540	src[14] = wk1r * (x0i - x0r);
541	src[15] = wk1r * (x0i + x0r);
542	k1 = 0;
543	for (j = 16; j < _length; j += 16) {
544	k1 += 2;
545	k2 = 2 * k1;
546	wk2r = _waveTable[k1];
547	wk2i = _waveTable[k1+1];
548	wk1r = _waveTable[k2];
549	wk1i = _waveTable[k2+1];
550	wk3r = wk1r - 2 * wk2i * wk1i;
551	wk3i = 2 * wk2i * wk1r - wk1i;
552	x0r = src[j] + src[j + 2];
553	x0i = src[j+1] + src[j + 3];
554	x1r = src[j] - src[j + 2];
555	x1i = src[j+1] - src[j + 3];
556	x2r = src[j + 4] + src[j + 6];
557	x2i = src[j + 5] + src[j + 7];
558	x3r = src[j + 4] - src[j + 6];
559	x3i = src[j + 5] - src[j + 7];
560	src[j] = x0r + x2r;
561	src[j+1] = x0i + x2i;
562	x0r -= x2r;
563	x0i -= x2i;
564	src[j + 4] = wk2r * x0r - wk2i * x0i;
565	src[j + 5] = wk2r * x0i + wk2i * x0r;
566	x0r = x1r - x3i;
567	x0i = x1i + x3r;
568	src[j + 2] = wk1r * x0r - wk1i * x0i;
569	src[j + 3] = wk1r * x0i + wk1i * x0r;
570	x0r = x1r + x3i;
571	x0i = x1i - x3r;
572	src[j + 6] = wk3r * x0r - wk3i * x0i;
573	src[j + 7] = wk3r * x0i + wk3i * x0r;
574	wk1r = _waveTable[k2 + 2];
575	wk1i = _waveTable[k2 + 3];
576	wk3r = wk1r - 2 * wk2r * wk1i;
577	wk3i = 2 * wk2r * wk1r - wk1i;
578	x0r = src[j + 8] + src[j+10];
579	x0i = src[j + 9] + src[j+11];
580	x1r = src[j + 8] - src[j+10];
581	x1i = src[j + 9] - src[j+11];
582	x2r = src[j+12] + src[j+14];
583	x2i = src[j+13] + src[j+15];
584	x3r = src[j+12] - src[j+14];
585	x3i = src[j+13] - src[j+15];
586	src[j + 8] = x0r + x2r;
587	src[j + 9] = x0i + x2i;
588	x0r -= x2r;
589	x0i -= x2i;
590	src[j+12] = -wk2i * x0r - wk2r * x0i;
591	src[j+13] = -wk2i * x0i + wk2r * x0r;
592	x0r = x1r - x3i;
593	x0i = x1i + x3r;
594	src[j+10] = wk1r * x0r - wk1i * x0i;
595	src[j+11] = wk1r * x0i + wk1i * x0r;
596	x0r = x1r + x3i;
597	x0i = x1i - x3r;
598	src[j+14] = wk3r * x0r - wk3i * x0i;
599	src[j+15] = wk3r * x0i + wk3i * x0r;
600	}
601	}
602
603
604	private function _cftmdl(src:Vector.<Number>, l:int) : void
605	{
606	var j:int, j1:int, j2:int, j3:int, k:int, k1:int, k2:int, m:int, m2:int,
607	wk1r:Number, wk2r:Number, wk3r:Number, x0r:Number, x1r:Number, x2r:Number, x3r:Number,
608	wk1i:Number, wk2i:Number, wk3i:Number, x0i:Number, x1i:Number, x2i:Number, x3i:Number;
609
610	m = l << 2;
611	for (j = 0; j < l; j += 2) {
612	j1 = j + l;
613	j2 = j1 + l;
614	j3 = j2 + l;
615	x0r = src[j] + src[j1];
616	x0i = src[j+1] + src[j1+1];
617	x1r = src[j] - src[j1];
618	x1i = src[j+1] - src[j1+1];
619	x2r = src[j2] + src[j3];
620	x2i = src[j2+1] + src[j3+1];
621	x3r = src[j2] - src[j3];
622	x3i = src[j2+1] - src[j3+1];
623	src[j] = x0r + x2r;
624	src[j+1] = x0i + x2i;
625	src[j2] = x0r - x2r;
626	src[j2+1] = x0i - x2i;
627	src[j1] = x1r - x3i;
628	src[j1+1] = x1i + x3r;
629	src[j3] = x1r + x3i;
630	src[j3+1] = x1i - x3r;
631	}
632	wk1r = _waveTable[2];
633	for (j = m; j < l + m; j += 2) {
634	j1 = j + l;
635	j2 = j1 + l;
636	j3 = j2 + l;
637	x0r = src[j] + src[j1];
638	x0i = src[j+1] + src[j1+1];
639	x1r = src[j] - src[j1];
640	x1i = src[j+1] - src[j1+1];
641	x2r = src[j2] + src[j3];
642	x2i = src[j2+1] + src[j3+1];
643	x3r = src[j2] - src[j3];
644	x3i = src[j2+1] - src[j3+1];
645	src[j] = x0r + x2r;
646	src[j+1] = x0i + x2i;
647	src[j2] = x2i - x0i;
648	src[j2+1] = x0r - x2r;
649	x0r = x1r - x3i;
650	x0i = x1i + x3r;
651	src[j1] = wk1r * (x0r - x0i);
652	src[j1+1] = wk1r * (x0r + x0i);
653	x0r = x3i + x1r;
654	x0i = x3r - x1i;
655	src[j3] = wk1r * (x0i - x0r);
656	src[j3+1] = wk1r * (x0i + x0r);
657	}
658	k1 = 0;
659	m2 = 2 * m;
660	for (k = m2; k < _length; k += m2) {
661	k1 += 2;
662	k2 = 2 * k1;
663	wk2r = _waveTable[k1];
664	wk2i = _waveTable[k1+1];
665	wk1r = _waveTable[k2];
666	wk1i = _waveTable[k2+1];
667	wk3r = wk1r - 2 * wk2i * wk1i;
668	wk3i = 2 * wk2i * wk1r - wk1i;
669	for (j = k; j < l + k; j += 2) {
670	j1 = j + l;
671	j2 = j1 + l;
672	j3 = j2 + l;
673	x0r = src[j] + src[j1];
674	x0i = src[j+1] + src[j1+1];
675	x1r = src[j] - src[j1];
676	x1i = src[j+1] - src[j1+1];
677	x2r = src[j2] + src[j3];
678	x2i = src[j2+1] + src[j3+1];
679	x3r = src[j2] - src[j3];
680	x3i = src[j2+1] - src[j3+1];
681	src[j] = x0r + x2r;
682	src[j+1] = x0i + x2i;
683	x0r -= x2r;
684	x0i -= x2i;
685	src[j2] = wk2r * x0r - wk2i * x0i;
686	src[j2+1] = wk2r * x0i + wk2i * x0r;
687	x0r = x1r - x3i;
688	x0i = x1i + x3r;
689	src[j1] = wk1r * x0r - wk1i * x0i;
690	src[j1+1] = wk1r * x0i + wk1i * x0r;
691	x0r = x1r + x3i;
692	x0i = x1i - x3r;
693	src[j3] = wk3r * x0r - wk3i * x0i;
694	src[j3+1] = wk3r * x0i + wk3i * x0r;
695	}
696	wk1r = _waveTable[k2 + 2];
697	wk1i = _waveTable[k2 + 3];
698	wk3r = wk1r - 2 * wk2r * wk1i;
699	wk3i = 2 * wk2r * wk1r - wk1i;
700	for (j = k + m; j < l + (k + m); j += 2) {
701	j1 = j + l;
702	j2 = j1 + l;
703	j3 = j2 + l;
704	x0r = src[j] + src[j1];
705	x0i = src[j+1] + src[j1+1];
706	x1r = src[j] - src[j1];
707	x1i = src[j+1] - src[j1+1];
708	x2r = src[j2] + src[j3];
709	x2i = src[j2+1] + src[j3+1];
710	x3r = src[j2] - src[j3];
711	x3i = src[j2+1] - src[j3+1];
712	src[j] = x0r + x2r;
713	src[j+1] = x0i + x2i;
714	x0r -= x2r;
715	x0i -= x2i;
716	src[j2] = -wk2i * x0r - wk2r * x0i;
717	src[j2+1] = -wk2i * x0i + wk2r * x0r;
718	x0r = x1r - x3i;
719	x0i = x1i + x3r;
720	src[j1] = wk1r * x0r - wk1i * x0i;
721	src[j1+1] = wk1r * x0i + wk1i * x0r;
722	x0r = x1r + x3i;
723	x0i = x1i - x3r;
724	src[j3] = wk3r * x0r - wk3i * x0i;
725	src[j3+1] = wk3r * x0i + wk3i * x0r;
726	}
727	}
728	}
729
730
731	private function _rftfsub(src:Vector.<Number>) : void
732	{
733	var j:int, k:int, kk:int, m:int,
734	wkr:Number, wki:Number, xr:Number, xi:Number, yr:Number, yi:Number;
735
736	m = _length >> 1;
737	kk = 0;
738	for (j = 2; j < m; j += 2) {
739	k = _length - j;
740	kk += 4;
741	wkr = 0.5 - _cosTable[_length - kk];
742	wki = _cosTable[kk];
743	xr = src[j] - src[k];
744	xi = src[j+1] + src[k+1];
745	yr = wkr * xr - wki * xi;
746	yi = wkr * xi + wki * xr;
747	src[j] -= yr;
748	src[j+1] -= yi;
749	src[k] += yr;
750	src[k+1] -= yi;
751	}
752	}
753
754
755	private function _rftbsub(src:Vector.<Number>) : void
756	{
757	var j:int, k:int, kk:int, m:int,
758	wkr:Number, wki:Number, xr:Number, xi:Number, yr:Number, yi:Number;
759
760	src[1] = -src[1];
761	m = _length >> 1;
762	kk = 0;
763	for (j = 2; j < m; j += 2) {
764	k = _length - j;
765	kk += 4;
766	wkr = 0.5 - _cosTable[_length - kk];
767	wki = _cosTable[kk];
768	xr = src[j] - src[k];
769	xi = src[j+1] + src[k+1];
770	yr = wkr * xr + wki * xi;
771	yi = wkr * xi - wki * xr;
772	src[j] -= yr;
773	src[j+1] = yi - src[j+1];
774	src[k] += yr;
775	src[k+1] = yi - src[k+1];
776	}
777	src[m+1] = -src[m+1];
778	}
779
780
781	private function _dctsub(src:Vector.<Number>) : void
782	{
783	var j:int, k:int, kk:int, m:int, wkr:Number, wki:Number, xr:Number;
784
785	m = _length >> 1;
786	kk = 0;
787	for (j = 1; j < m; j++) {
788	k = _length - j;
789	kk += 1;
790	wkr = _cosTable[kk] - _cosTable[_length - kk];
791	wki = _cosTable[kk] + _cosTable[_length - kk];
792	xr = wki * src[j] - wkr * src[k];
793	src[j] = wkr * src[j] + wki * src[k];
794	src[k] = xr;
795	}
796	src[m] *= 0.7071067811865476; // cos(pi/4)
797	}
798
799
800	private function _dstsub(src:Vector.<Number>) : void
801	{
802	var j:int, k:int, kk:int, m:int, wkr:Number, wki:Number, xr:Number;
803
804	m = _length >> 1;
805	kk = 0;
806	for (j = 1; j < m; j++) {
807	k = _length - j;
808	kk += 1;
809	wkr = _cosTable[kk] - _cosTable[_length - kk];
810	wki = _cosTable[kk] + _cosTable[_length - kk];
811	xr = wki * src[k] - wkr * src[j];
812	src[k] = wkr * src[k] + wki * src[j];
813	src[j] = xr;
814	}
815	src[m] *= 0.7071067811865476; // cos(pi/4)
816	}
817	}
818	}
819
820
821
822
823
824
825
826
827
828

Note: リポジトリブラウザについてのヘルプは TracBrowser を参照してください。

異なるフォーマットでダウンロード:

Original Format

Spark project

root/as3/SiOPM/trunk/src/org/si/utils/FFT_original.as

異なるフォーマットでダウンロード: