1 /*
2 * real.c
3 *
4 * Created on: Apr 20, 2013
5 * Author: Rafat Hussain
6 */
7 #include <stdio.h>
8 #include "real.h"
9
fft_real_init(int N,int sgn)10 fft_real_object fft_real_init(int N, int sgn) {
11 fft_real_object obj = NULL;
12 fft_type theta;
13 int k;
14
15 obj = (fft_real_object) malloc (sizeof(struct fft_real_set) + sizeof(fft_data)* (N/2));
16
17 obj->cobj = fft_init(N/2,sgn);
18
19 for (k = 0; k < N/2;++k) {
20 theta = PI2*k/N;
21 obj->twiddle2[k].re = cos(theta);
22 obj->twiddle2[k].im = sin(theta);
23
24 }
25
26
27 return obj;
28
29
30 }
31
fft_r2c_exec(fft_real_object obj,fft_type * inp,fft_data * oup)32 void fft_r2c_exec(fft_real_object obj,fft_type *inp,fft_data *oup) {
33 fft_data* cinp;
34 fft_data* coup;
35 int i,N2,N;
36 fft_type temp1,temp2;
37 N2 = obj->cobj->N;
38 N = N2*2;
39
40 cinp = (fft_data*) malloc (sizeof(fft_data) * N2);
41 coup = (fft_data*) malloc (sizeof(fft_data) * N2);
42
43 for (i = 0; i < N2; ++i) {
44 cinp[i].re = inp[2*i];
45 cinp[i].im = inp[2*i+1];
46 }
47
48 fft_exec(obj->cobj,cinp,coup);
49
50 oup[0].re = coup[0].re + coup[0].im;
51 oup[0].im = 0.0;
52
53 for (i = 1; i < N2; ++i) {
54 temp1 = coup[i].im + coup[N2-i].im ;
55 temp2 = coup[N2-i].re - coup[i].re ;
56 oup[i].re = (coup[i].re + coup[N2-i].re + (temp1 * obj->twiddle2[i].re) + (temp2 * obj->twiddle2[i].im)) / 2.0;
57 oup[i].im = (coup[i].im - coup[N2-i].im + (temp2 * obj->twiddle2[i].re) - (temp1 * obj->twiddle2[i].im)) / 2.0;
58 }
59
60
61
62 oup[N2].re = coup[0].re - coup[0].im;
63 oup[N2].im = 0.0;
64
65 for (i = 1; i < N2;++i) {
66 oup[N-i].re = oup[i].re;
67 oup[N-i].im = -oup[i].im;
68 }
69
70
71 free(cinp);
72 free(coup);
73
74 }
75
fft_c2r_exec(fft_real_object obj,fft_data * inp,fft_type * oup)76 void fft_c2r_exec(fft_real_object obj,fft_data *inp,fft_type *oup) {
77
78 fft_data* cinp;
79 fft_data* coup;
80 int i,N2;
81 fft_type temp1,temp2;
82 N2 = obj->cobj->N;
83
84 cinp = (fft_data*) malloc (sizeof(fft_data) * N2);
85 coup = (fft_data*) malloc (sizeof(fft_data) * N2);
86
87 for (i = 0; i < N2; ++i) {
88 temp1 = -inp[i].im - inp[N2-i].im ;
89 temp2 = -inp[N2-i].re + inp[i].re ;
90 cinp[i].re = inp[i].re + inp[N2-i].re + (temp1 * obj->twiddle2[i].re) - (temp2 * obj->twiddle2[i].im);
91 cinp[i].im = inp[i].im - inp[N2-i].im + (temp2 * obj->twiddle2[i].re) + (temp1 * obj->twiddle2[i].im);
92 }
93
94 fft_exec(obj->cobj,cinp,coup);
95 for (i = 0; i < N2; ++i) {
96 oup[2*i] = coup[i].re;
97 oup[2*i+1] = coup[i].im;
98 }
99 free(cinp);
100 free(coup);
101
102
103 }
104
free_real_fft(fft_real_object object)105 void free_real_fft(fft_real_object object) {
106 free_fft(object->cobj);
107 free(object);
108 }
109
110