DFT/dft.h at main · code4fukui/DFT · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
#include <math.h>

#ifndef __CUDACC__
#define __device__
#endif

void dft(float* src, size_t len, float* re, float* im) {
  for (int i = 0; i < len; i++) {
    float re1 = 0.0f;
    float im1 = 0.0f;
    for (int j = 0; j < len; j++) {
      float th = (2.0f * (float)M_PI) * i * j / len;
      float d = src[j];
      re1 += d * cosf(th);
      im1 -= d * sinf(th);
    }
    re[i] = re1;
    im[i] = im1;
  }
  /*
  for (int i = 0; i < len; i++) {
    float re1 = 0.0f;
    float im1 = 0.0f;
    float dth = (2.0f * (float)M_PI) * i / len;
    float th = 0.0f;
    for (int j = 0; j < len; j++) {
      float d = src[j];
      re1 += d * cosf(th);
      im1 -= d * sinf(th);
      th += dth;
    }
    re[i] = re1;
    im[i] = im1;
  }
  */
}

void idft(float* re, float* im, size_t len, float* dst) {
  for (int i = 0; i < len; i++) {
    float re1 = 0.0f;
    for (int j = 0; j < len; j++) {
      float th = (2.0f * (float)M_PI) * i * j / len;
      re1 += re[j] * cosf(th) - im[j] * sinf(th);
      //im1 += re[j] * Math.sin(th) + im[j] * Math.cos(th);
    }
    dst[i] = re1 / len;
    //im2[i] = im1 / len;
  }

  /*
  for (int i = 0; i < len; i++) {
    float re1 = 0.0f;
    float dth = (2.0f * (float)M_PI) * i / len;
    float th = 0.0f;
    for (int j = 0; j < len; j++) {
      re1 += re[j] * cosf(th) - im[j] * sinf(th);
      //im1 += re[j] * Math.sin(th) + im[j] * Math.cos(th);
      th += dth;
    }
    dst[i] = re1 / len;
    //im2[i] = im1 / len;
  }
  */

  /*
  for (int i = 0; i < len; i++) {
    float re1 = 0.0f;
    float th = 0.0f;
    float dth = (2.0f * (float)M_PI) * i / len;
    for (int j = 0; j < len; j++) {
      re1 += re[j] * cosf(th) - im[j] * sinf(th);
      th += dth;
    }
    dst[i] = re1 / len;
  }
  */
}

__device__ void dft1(float* src, size_t len, float* re, float* im, int i) {
  float re1 = 0.0f;
  float im1 = 0.0f;
  for (int j = 0; j < len; j++) {
    float th = (2.0f * (float)M_PI) * i * j/ len;
    float d = src[j];
    re1 += d * cosf(th);
    im1 -= d * sinf(th);
  }
  re[i] = re1;
  im[i] = im1;
  /*
  float re1 = 0.0f;
  float im1 = 0.0f;
  float dth = (2.0f * (float)M_PI) * i / len;
  float th = 0.0f;
  for (int j = 0; j < len; j++) {
    float d = src[j];
    re1 += d * cosf(th);
    im1 -= d * sinf(th);
    th += dth;
  }
  re[i] = re1;
  im[i] = im1;
  */
}