谁有图像 FFT 正反变换的源码

zhujiang73

谁有图像 FFT 正反变换的源码，最近在恶补图像处理。  {:3_189:}

mfzz1134

直接用FFTW就可以了，但要注意处理的时候如果有卷积可能会有移位的问题，可能需要扩充和补零。

cjaizss

正反几乎一样的,就一个系数差.
2^n*2^n的fft
先每个2*2变换,然后两两个两个根据公式合并

zhujiang73

正反几乎一样的,就一个系数差.
2^n*2^n的fft
先每个2*2变换,然后两两个两个根据公式合并
cjaizss 发表于 2010-12-28 14:11

    是不是要先算一维的FFT,还是有更好的算法。

cjaizss

是不是要先算一维的FFT,还是有更好的算法。
zhujiang73 发表于 2010-12-28 14:40

   还是二维的,你可以自己去查一下fft的递归公式

net_robber

楼上果然钟爱数学

id_for_fun

有库可用

zhujiang73

还是二维的,你可以自己去查一下fft的递归公式
cjaizss 发表于 2010-12-28 14:54

   是这样的？  google  到这个：  http://blog.163.com/ld081055@126 ... 691520105144852433/

1. #include <stdio.h>
   
2. #include <stdlib.h>
   
3. #include <math.h>
   
4. 
   
5. #define intsize sizeof(int)
   
6. #define complexsize sizeof(complex)
   
7. #define PI 3.1415926
   
8. 
   
9. int *a,*b;
   
10. int nLen,init_nLen,mLen,init_mLen,N,M;
    
11. FILE *dataFile;
    
12. 
    
13. typedef struct{
    
14.     float real;
    
15.     float image;
    
16. }complex;
    
17. 
    
18. complex *A,*A_In,*W;
    
19. 
    
20. complex Add(complex, complex);
    
21. complex Sub(complex, complex);
    
22. complex Mul(complex, complex);
    
23. int calculate_M(int);
    
24. void reverse(int,int);
    
25. void readData();
    
26. void fft(int,int);
    
27. void Ifft();
    
28. void printResult_fft();
    
29. void printResult_Ifft();
    
30. 
    
31. int main()
    
32. {
    
33.     int i,j;
    
34.     readData();
    
35.     A = (complex *)malloc(complexsize*nLen);
    
36.     reverse(nLen,N);
    
37.     for(i=0; i<mLen; i++)
    
38.     {
    
39.         for(j=0; j<nLen; j++)
    
40.         {
    
41.             A[j].real = A_In[i*nLen+b[j]].real;
    
42.             A[j].image = A_In[i*nLen+b[j]].image;
    
43.         }
    
44.       
    
45.         fft(nLen,N);
    
46.         for(j=0; j<nLen; j++)
    
47.         {
    
48.             A_In[i*nLen+j].real = A[j].real;
    
49.             A_In[i*nLen+j].image = A[j].image;
    
50.         }
    
51.     }
    
52. 
    
53.     free(a);
    
54.     free(b);
    
55.     free(A);
    
56. 
    
57.     A = (complex *)malloc(complexsize*mLen);
    
58.     reverse(mLen,M);
    
59.     for(i=0; i<nLen; i++)
    
60.     {
    
61.         for(j=0; j<mLen; j++)
    
62.         {
    
63.             A[j].real = A_In[b[j]*nLen+i].real;
    
64.             A[j].image = A_In[b[j]*nLen+i].image;
    
65.         }
    
66. 
    
67.         fft(mLen,M);
    
68.         for(j=0; j<mLen; j++)
    
69.         {
    
70.             A_In[j*nLen+i].real = A[j].real;
    
71.             A_In[j*nLen+i].image = A[j].image;
    
72.         }
    
73.     }
    
74.     free(A);
    
75.     printResult_fft();
    
76.     Ifft();
    
77.     printResult_Ifft();
    
78.     return 0;
    
79. }
    
80. 
    
81. void readData()
    
82. {
    
83.      int i,j;
    
84. 
    
85.      dataFile = fopen("dataIn.txt","r");
    
86.      fscanf(dataFile,"%d %d",&init_mLen,&init_nLen);
    
87.      M = calculate_M(init_mLen);
    
88.      N = calculate_M(init_nLen);
    
89.      nLen = (int)pow(2,N);
    
90.      mLen = (int)pow(2,M);
    
91.      A_In = (complex *)malloc(complexsize*nLen*mLen);
    
92. 
    
93.      for(i=0; i<init_mLen; i++)
    
94.      {
    
95.          for(j=0; j<init_nLen; j++)
    
96.          {
    
97.              fscanf(dataFile,"%f",&A_In[i*nLen+j].real);
    
98.              A_In[i*nLen+j].image = 0.0;
    
99.          }
    
100.      }
     
101.      fclose(dataFile);
     
102. 
     
103.      for(i=0; i<mLen; i++)
     
104.      {
     
105.          for(j=init_nLen; j<nLen; j++)
     
106.          {
     
107.              A_In[i*nLen+j].real = 0.0;
     
108.              A_In[i*nLen+j].image = 0.0;
     
109.          }
     
110.      }
     
111. 
     
112.      for(i=init_mLen; i<mLen; i++)
     
113.      {
     
114.          for(j=0; j<init_nLen; j++)
     
115.          {
     
116.              A_In[i*nLen+j].real = 0.0;
     
117.              A_In[i*nLen+j].image = 0.0;
     
118.          }
     
119.      }
     
120. 
     
121.      printf("Reading initial datas:\n");
     
122.      for(i=0; i<init_mLen; i++)
     
123.      {
     
124.          for(j=0; j<init_nLen; j++)
     
125.          {
     
126.              if(A_In[i*nLen+j].image < 0)
     
127.              {
     
128.                  printf("%f%fi\t",A_In[i*nLen+j].real,A_In[i*nLen+j].image);
     
129.              }
     
130.              else
     
131.              {
     
132.                  printf("%f+%fi\t",A_In[i*nLen+j].real,A_In[i*nLen+j].image);
     
133.              }
     
134.          }
     
135.          printf("\n");
     
136.      }
     
137. 
     
138.      printf("\n");
     
139. 
     
140.      printf("Reading formal datas:\n");
     
141.      for(i=0; i<mLen; i++)
     
142.      {
     
143.          for(j=0; j<nLen; j++)
     
144.          {
     
145.              if(A_In[i*nLen+j].image < 0)
     
146.              {
     
147.                  printf("%f%fi\t",A_In[i*nLen+j].real,A_In[i*nLen+j].image);
     
148.              }
     
149.              else
     
150.              {
     
151.                  printf("%f+%fi\t",A_In[i*nLen+j].real,A_In[i*nLen+j].image);
     
152.              }
     
153.          }
     
154.          printf("\n");
     
155.      }
     
156. }
     
157. 
     
158. 
     
159. void fft(int fft_nLen, int fft_M)
     
160. {
     
161.      int i;
     
162.      int lev,dist,p,t;
     
163.      complex B;
     
164. 
     
165.      W = (complex *)malloc(complexsize*fft_nLen/2);
     
166. 
     
167.      for(lev=1; lev<=fft_M; lev++)
     
168.      {
     
169.          dist = (int)pow(2,lev-1);
     
170.          for(t=0; t<dist; t++)
     
171.          {
     
172.              p = t*(int)pow(2,fft_M-lev);
     
173.              W[p].real = (float)cos(2*PI*p/fft_nLen);
     
174.              W[p].image = (float)(-1*sin(2*PI*p/fft_nLen));
     
175.              for(i=t; i<fft_nLen; i=i+(int)pow(2,lev))
     
176.              {
     
177.                  B = Add(A[i],Mul(A[i+dist],W[p]));
     
178.                  A[i+dist] = Sub(A[i],Mul(A[i+dist],W[p]));
     
179.                  A[i].real = B.real;
     
180.                  A[i].image = B.image;
     
181.              }
     
182.          }
     
183.      }
     
184. 
     
185.      free(W);
     
186. }
     
187. 
     
188. 
     
189. void printResult_fft()
     
190. {
     
191.      int i,j;
     
192. 
     
193.      printf("Output FFT results:\n");
     
194.      for(i=0; i<mLen; i++)
     
195.      {
     
196.          for(j=0; j<nLen; j++)
     
197.          {
     
198.              if(A_In[i*nLen+j].image < 0)
     
199.              {
     
200.                  printf("%f%fi\t",A_In[i*nLen+j].real,A_In[i*nLen+j].image);
     
201.              }
     
202.              else
     
203.              {
     
204.                  printf("%f+%fi\t",A_In[i*nLen+j].real,A_In[i*nLen+j].image);
     
205.              }
     
206.          }
     
207.          printf("\n");
     
208.      }
     
209. }
     
210. 
     
211. void printResult_Ifft()
     
212. {
     
213.      int i,j;
     
214. 
     
215.      printf("Output IFFT results:\n");
     
216.      for(i=0; i<mLen; i++)
     
217.      {
     
218.          for(j=0; j<nLen; j++)
     
219.          {
     
220.              if(A_In[i*nLen+j].image < 0)
     
221.              {
     
222.                  printf("%f%fi\t",A_In[i*nLen+j].real,A_In[i*nLen+j].image);
     
223.              }
     
224.              else
     
225.              {
     
226.                  printf("%f+%fi\t",A_In[i*nLen+j].real,A_In[i*nLen+j].image);
     
227.              }
     
228.          }
     
229.          printf("\n");
     
230.      }
     
231. 
     
232.      free(A_In);
     
233. }
     
234. 
     
235. int calculate_M(int len)
     
236. {
     
237.     int i;
     
238.     int k;
     
239. 
     
240.     i = 0;
     
241.     k = 1;
     
242.     while(k < len)
     
243.     {
     
244.         k = k*2;
     
245.         i++;
     
246.     }
     
247. 
     
248.     return i;
     
249. }
     
250. 
     
251. void reverse(int len, int M)
     
252. {
     
253.     int i,j;
     
254. 
     
255.     a = (int *)malloc(intsize*M);
     
256.     b = (int *)malloc(intsize*len);
     
257. 
     
258.     for(i=0; i<M; i++)
     
259.     {
     
260.         a[i] = 0;
     
261.     }
     
262. 
     
263.     b[0] = 0;
     
264.     for(i=1; i<len; i++)
     
265.     {
     
266.         j = 0;
     
267.         while(a[j] != 0)
     
268.         {
     
269.             a[j] = 0;
     
270.             j++;
     
271.         }
     
272. 
     
273.         a[j] = 1;
     
274.         b[i] = 0;
     
275.         for(j=0; j<M; j++)
     
276.         {
     
277.             b[i] = b[i]+a[j]*(int)pow(2,M-1-j);
     
278.         }
     
279.     }
     
280. }
     
281. 
     
282. complex Add(complex c1, complex c2)
     
283. {
     
284.     complex c;
     
285.     c.real = c1.real+c2.real;
     
286.     c.image = c1.image+c2.image;
     
287.     return c;
     
288. }
     
289. 
     
290. complex Sub(complex c1, complex c2)
     
291. {
     
292.     complex c;
     
293.     c.real = c1.real-c2.real;
     
294.     c.image = c1.image-c2.image;
     
295.     return c;
     
296. }
     
297. 
     
298. complex Mul(complex c1, complex c2)
     
299. {
     
300.     complex c;
     
301.     c.real = c1.real*c2.real-c1.image*c2.image;
     
302.     c.image = c1.real*c2.image+c2.real*c1.image;
     
303.     return c;
     
304. }
     
305. 
     
306. void Ifft()
     
307. {
     
308.     int i,j;
     
309. 
     
310.     for(i=0; i<mLen; i++)
     
311.     {
     
312.         for(j=0; j<nLen; j++)
     
313.         {
     
314.             A_In[i*nLen+j].image = -A_In[i*nLen+j].image;
     
315.         }
     
316.     }
     
317. 
     
318.     A = (complex *)malloc(complexsize*nLen);
     
319.     reverse(nLen,N);
     
320.     for(i=0; i<mLen; i++)
     
321.     {
     
322.         for(j=0; j<nLen; j++)
     
323.         {
     
324.             A[j].real = A_In[i*nLen+b[j]].real;
     
325.             A[j].image = A_In[i*nLen+b[j]].image;  
     
326.         }
     
327. 
     
328.         fft(nLen,N);
     
329.         for(j=0; j<nLen; j++)
     
330.         {   
     
331.             A_In[i*nLen+j].real = A[j].real/nLen;
     
332.             A_In[i*nLen+j].image = A[j].image/nLen;
     
333.         }
     
334.     }
     
335.     free(A);
     
336.     free(a);
     
337.     free(b);
     
338. 
     
339.     A = (complex *)malloc(complexsize*mLen);
     
340.     reverse(mLen,M);
     
341.     for(i=0; i<nLen; i++)
     
342.     {
     
343.         for(j=0; j<mLen; j++)
     
344.         {
     
345.             A[j].real = A_In[b[j]*nLen+i].real;
     
346.             A[j].image = A_In[b[j]*nLen+i].image;
     
347.         }
     
348. 
     
349.         fft(mLen,M);
     
350.         for(j=0; j<mLen; j++)
     
351.         {
     
352.             A_In[j*nLen+i].real = A[j].real/mLen;
     
353.             A_In[j*nLen+i].image = A[j].image/mLen;
     
354.         }
     
355.     }
     
356.     free(A);
     
357.     free(a);
     
358.     free(b);
     
359. }

复制代码

zhujiang73

直接用FFTW就可以了，但要注意处理的时候如果有卷积可能会有移位的问题，可能需要扩充和补零。
mfzz1134 发表于 2010-12-29 13:45

    是要算卷积用，我想做个卷积的 GPU 加速，想找个 FFT 的源码参考参考...。  {:3_189:}