API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Dec 4th, 2014
100
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 7.08 KB | None | 0 0
raw download clone embed print report
  1. #include <avr/pgmspace.h>
  2. #include "fix_fft.h"
  3. #include <Arduino.h>
  4.  
  5. /* fix_fft.c - Fixed-point in-place Fast Fourier Transform */
  6. /*
  7. All data are fixed-point short integers, in which -32768
  8. to +32768 represent -1.0 to +1.0 respectively. Integer
  9. arithmetic is used for speed, instead of the more natural
  10. floating-point.
  11.  
  12. For the forward FFT (time -> freq), fixed scaling is
  13. performed to prevent arithmetic overflow, and to map a 0dB
  14. sine/cosine wave (i.e. amplitude = 32767) to two -6dB freq
  15. coefficients. The return value is always 0.
  16.  
  17. For the inverse FFT (freq -> time), fixed scaling cannot be
  18. done, as two 0dB coefficients would sum to a peak amplitude
  19. of 64K, overflowing the 32k range of the fixed-point integers.
  20. Thus, the fix_fft() routine performs variable scaling, and
  21. returns a value which is the number of bits LEFT by which
  22. the output must be shifted to get the actual amplitude
  23. (i.e. if fix_fft() returns 3, each value of fr[] and fi[]
  24. must be multiplied by 8 (2**3) for proper scaling.
  25. Clearly, this cannot be done within fixed-point short
  26. integers. In practice, if the result is to be used as a
  27. filter, the scale_shift can usually be ignored, as the
  28. result will be approximately correctly normalized as is.
  29.  
  30. Written by: Tom Roberts 11/8/89
  31. Made portable: Malcolm Slaney 12/15/94 malcolm@interval.com
  32. Enhanced: Dimitrios P. Bouras 14 Jun 2006 dbouras@ieee.org
  33. Modified for 8bit values David Keller 10.10.2010
  34. */
  35.  
  36.  
  37. #define N_WAVE 256 /* full length of Sinewave[] */
  38. #define LOG2_N_WAVE 8 /* log2(N_WAVE) */
  39.  
  40.  
  41.  
  42.  
  43. /*
  44. Since we only use 3/4 of N_WAVE, we define only
  45. this many samples, in order to conserve data space.
  46. */
  47.  
  48.  
  49.  
  50. const prog_int8_t Sinewave[N_WAVE-N_WAVE/4] PROGMEM = {
  51. 0, 3, 6, 9, 12, 15, 18, 21,
  52. 24, 28, 31, 34, 37, 40, 43, 46,
  53. 48, 51, 54, 57, 60, 63, 65, 68,
  54. 71, 73, 76, 78, 81, 83, 85, 88,
  55. 90, 92, 94, 96, 98, 100, 102, 104,
  56. 106, 108, 109, 111, 112, 114, 115, 117,
  57. 118, 119, 120, 121, 122, 123, 124, 124,
  58. 125, 126, 126, 127, 127, 127, 127, 127,
  59.  
  60. 127, 127, 127, 127, 127, 127, 126, 126,
  61. 125, 124, 124, 123, 122, 121, 120, 119,
  62. 118, 117, 115, 114, 112, 111, 109, 108,
  63. 106, 104, 102, 100, 98, 96, 94, 92,
  64. 90, 88, 85, 83, 81, 78, 76, 73,
  65. 71, 68, 65, 63, 60, 57, 54, 51,
  66. 48, 46, 43, 40, 37, 34, 31, 28,
  67. 24, 21, 18, 15, 12, 9, 6, 3,
  68.  
  69. 0, -3, -6, -9, -12, -15, -18, -21,
  70. -24, -28, -31, -34, -37, -40, -43, -46,
  71. -48, -51, -54, -57, -60, -63, -65, -68,
  72. -71, -73, -76, -78, -81, -83, -85, -88,
  73. -90, -92, -94, -96, -98, -100, -102, -104,
  74. -106, -108, -109, -111, -112, -114, -115, -117,
  75. -118, -119, -120, -121, -122, -123, -124, -124,
  76. -125, -126, -126, -127, -127, -127, -127, -127,
  77.  
  78. /*-127, -127, -127, -127, -127, -127, -126, -126,
  79. -125, -124, -124, -123, -122, -121, -120, -119,
  80. -118, -117, -115, -114, -112, -111, -109, -108,
  81. -106, -104, -102, -100, -98, -96, -94, -92,
  82. -90, -88, -85, -83, -81, -78, -76, -73,
  83. -71, -68, -65, -63, -60, -57, -54, -51,
  84. -48, -46, -43, -40, -37, -34, -31, -28,
  85. -24, -21, -18, -15, -12, -9, -6, -3, */
  86. };
  87.  
  88.  
  89.  
  90.  
  91.  
  92.  
  93. /*
  94. FIX_MPY() - fixed-point multiplication & scaling.
  95. Substitute inline assembly for hardware-specific
  96. optimization suited to a particluar DSP processor.
  97. Scaling ensures that result remains 16-bit.
  98. */
  99. inline char FIX_MPY(char a, char b)
  100. {
  101.  
  102. //Serial.println(a);
  103. //Serial.println(b);
  104.  
  105.  
  106. /* shift right one less bit (i.e. 15-1) */
  107. int c = ((int)a * (int)b) >> 6;
  108. /* last bit shifted out = rounding-bit */
  109. b = c & 0x01;
  110. /* last shift + rounding bit */
  111. a = (c >> 1) + b;
  112.  
  113. /*
  114. Serial.println(Sinewave[3]);
  115. Serial.println(c);
  116. Serial.println(a);
  117. while(1);*/
  118.  
  119. return a;
  120. }
  121.  
  122. /*
  123. fix_fft() - perform forward/inverse fast Fourier transform.
  124. fr[n],fi[n] are real and imaginary arrays, both INPUT AND
  125. RESULT (in-place FFT), with 0 <= n < 2**m; set inverse to
  126. 0 for forward transform (FFT), or 1 for iFFT.
  127. */
  128. int fix_fft(char fr[], char fi[], int m, int inverse)
  129. {
  130. int mr, nn, i, j, l, k, istep, n, scale, shift;
  131. char qr, qi, tr, ti, wr, wi;
  132.  
  133. n = 1 << m;
  134.  
  135. /* max FFT size = N_WAVE */
  136. if (n > N_WAVE)
  137. return -1;
  138.  
  139. mr = 0;
  140. nn = n - 1;
  141. scale = 0;
  142.  
  143. /* decimation in time - re-order data */
  144. for (m=1; m<=nn; ++m) {
  145. l = n;
  146. do {
  147. l >>= 1;
  148. } while (mr+l > nn);
  149. mr = (mr & (l-1)) + l;
  150.  
  151. if (mr <= m)
  152. continue;
  153. tr = fr[m];
  154. fr[m] = fr[mr];
  155. fr[mr] = tr;
  156. ti = fi[m];
  157. fi[m] = fi[mr];
  158. fi[mr] = ti;
  159. }
  160.  
  161. l = 1;
  162. k = LOG2_N_WAVE-1;
  163. while (l < n) {
  164. if (inverse) {
  165. /* variable scaling, depending upon data */
  166. shift = 0;
  167. for (i=0; i<n; ++i) {
  168. j = fr[i];
  169. if (j < 0)
  170. j = -j;
  171. m = fi[i];
  172. if (m < 0)
  173. m = -m;
  174. if (j > 16383 || m > 16383) {
  175. shift = 1;
  176. break;
  177. }
  178. }
  179. if (shift)
  180. ++scale;
  181. } else {
  182. /*
  183. fixed scaling, for proper normalization --
  184. there will be log2(n) passes, so this results
  185. in an overall factor of 1/n, distributed to
  186. maximize arithmetic accuracy.
  187. */
  188. shift = 1;
  189. }
  190. /*
  191. it may not be obvious, but the shift will be
  192. performed on each data point exactly once,
  193. during this pass.
  194. */
  195. istep = l << 1;
  196. for (m=0; m<l; ++m) {
  197. j = m << k;
  198. /* 0 <= j < N_WAVE/2 */
  199. wr = pgm_read_word_near(Sinewave + j+N_WAVE/4);
  200.  
  201. /*Serial.println("asdfasdf");
  202. Serial.println(wr);
  203. Serial.println(j+N_WAVE/4);
  204. Serial.println(Sinewave[256]);
  205.  
  206. Serial.println("");*/
  207.  
  208.  
  209. wi = -pgm_read_word_near(Sinewave + j);
  210. if (inverse)
  211. wi = -wi;
  212. if (shift) {
  213. wr >>= 1;
  214. wi >>= 1;
  215. }
  216. for (i=m; i<n; i+=istep) {
  217. j = i + l;
  218. tr = FIX_MPY(wr,fr[j]) - FIX_MPY(wi,fi[j]);
  219. ti = FIX_MPY(wr,fi[j]) + FIX_MPY(wi,fr[j]);
  220. qr = fr[i];
  221. qi = fi[i];
  222. if (shift) {
  223. qr >>= 1;
  224. qi >>= 1;
  225. }
  226. fr[j] = qr - tr;
  227. fi[j] = qi - ti;
  228. fr[i] = qr + tr;
  229. fi[i] = qi + ti;
  230. }
  231. }
  232. --k;
  233. l = istep;
  234. }
  235. return scale;
  236. }
  237.  
  238. /*
  239. fix_fftr() - forward/inverse FFT on array of real numbers.
  240. Real FFT/iFFT using half-size complex FFT by distributing
  241. even/odd samples into real/imaginary arrays respectively.
  242. In order to save data space (i.e. to avoid two arrays, one
  243. for real, one for imaginary samples), we proceed in the
  244. following two steps: a) samples are rearranged in the real
  245. array so that all even samples are in places 0-(N/2-1) and
  246. all imaginary samples in places (N/2)-(N-1), and b) fix_fft
  247. is called with fr and fi pointing to index 0 and index N/2
  248. respectively in the original array. The above guarantees
  249. that fix_fft "sees" consecutive real samples as alternating
  250. real and imaginary samples in the complex array.
  251. */
  252. int fix_fftr(char f[], int m, int inverse)
  253. {
  254. int i, N = 1<<(m-1), scale = 0;
  255. char tt, *fr=f, *fi=&f[N];
  256.  
  257. if (inverse)
  258. scale = fix_fft(fi, fr, m-1, inverse);
  259. for (i=1; i<N; i+=2) {
  260. tt = f[N+i-1];
  261. f[N+i-1] = f[i];
  262. f[i] = tt;
  263. }
  264. if (! inverse)
  265. scale = fix_fft(fi, fr, m-1, inverse);
  266. return scale;
  267. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!