----------------------------------------------------------------------------- name: rfft() desc: real value fft these routines from the CARL software, spect.c check out the CARL CMusic distribution for more source code if forward is true, rfft replaces 2*N real data points in x with N complex values representing the positive frequency half of their Fourier spectrum, with x[1] replaced with the r
| 113 | // |
| 114 | //----------------------------------------------------------------------------- |
| 115 | void rfft(float * x, long N, unsigned int forward) |
| 116 | { |
| 117 | static int first = 1; |
| 118 | float c1, c2, h1r, h1i, h2r, h2i, wr, wi, wpr, wpi, temp, theta; |
| 119 | float xr, xi; |
| 120 | long i, i1, i2, i3, i4, N2p1; |
| 121 | |
| 122 | if (first) |
| 123 | { |
| 124 | PI = (float)(4.*atan(1.)); |
| 125 | TWOPI = (float)(8.*atan(1.)); |
| 126 | first = 0; |
| 127 | } |
| 128 | |
| 129 | theta = PI / N; |
| 130 | wr = 1.; |
| 131 | wi = 0.; |
| 132 | c1 = 0.5; |
| 133 | |
| 134 | if (forward) |
| 135 | { |
| 136 | c2 = -0.5; |
| 137 | cfft(x, N, forward); |
| 138 | xr = x[0]; |
| 139 | xi = x[1]; |
| 140 | } |
| 141 | else |
| 142 | { |
| 143 | c2 = 0.5; |
| 144 | theta = -theta; |
| 145 | xr = x[1]; |
| 146 | xi = 0.; |
| 147 | x[1] = 0.; |
| 148 | } |
| 149 | |
| 150 | wpr = (float)(-2.*pow(sin(0.5*theta), 2.)); |
| 151 | wpi = (float)sin(theta); |
| 152 | N2p1 = (N << 1) + 1; |
| 153 | |
| 154 | for (i = 0; i <= N >> 1; i++) |
| 155 | { |
| 156 | i1 = i << 1; |
| 157 | i2 = i1 + 1; |
| 158 | i3 = N2p1 - i2; |
| 159 | i4 = i3 + 1; |
| 160 | if (i == 0) |
| 161 | { |
| 162 | h1r = c1*(x[i1] + xr); |
| 163 | h1i = c1*(x[i2] - xi); |
| 164 | h2r = -c2*(x[i2] + xi); |
| 165 | h2i = c2*(x[i1] - xr); |
| 166 | x[i1] = h1r + wr*h2r - wi*h2i; |
| 167 | x[i2] = h1i + wr*h2i + wi*h2r; |
| 168 | xr = h1r - wr*h2r + wi*h2i; |
| 169 | xi = -h1i + wr*h2i + wi*h2r; |
| 170 | } |
| 171 | else |
| 172 | { |