长度随机的信号FFT如何实现

lebronze

想用C语言实现FFT、IFFT，但是遇到的问题是：如果信号的长度不是2的整数幂，那应该如何处理？
尝试了下面两种方法：
（1）截断原信号到2的整数幂。
（2）在原信号的后面补零到2的整数幂。
但发现结果都不正确。
这个应该有成熟的处理方法吧，希望懂的大神指导指导，多谢

lebronze

多谢各位回复，目前的解决方法是：对原始信号末尾补零到2的整数幂，在进行FFT。
时域补零，相当于频域插值，不会改变频谱的包络。

impulse

直接调用FFTW库-几个公认的比较好的、效率比较高的FFT算法之一
http://www.fftw.org/

lebronze

impulse 发表于 2016-7-5 22:34
直接调用FFTW库-几个公认的比较好的、效率比较高的FFT算法之一
http://www.fftw.org/

多谢回复，不过我这边是要搞嵌入式开发，FFT，IFFT只能自己写，不能调用库

hcharlie

从编程角度看，FFT，IFFT 的精髓就是必须是2的整数幂，加零等方法都是凑数的近似方法。

eastar

参考fft matlab代码改成c行吗

impulse

eastar 发表于 2016-7-6 08:18
参考fft matlab代码改成c行吗

matlab fft用的是FFTW

TestGuru

C语言算法可参考： http://forum.vibunion.com/thread-136567-1-5.html

ZH----过客

这是我找到的C++的代码，不知道有没有用  你看一下

1. 时域抽取法FFT的C++实现。
   
2.      下面是算法的流程图
   
3. 
   
4.      倒序的流程图
   
5. 
   
6. C++实现代码：
   
7. 1. FFT.h
   
8. #pragma once
   
9. #ifndef FFT_H
   
10. #define FFT_H
    
11. #include <cstring>
    
12. #include <cmath>
    
13. using namespace std;
    
14. #define pi 3.141592
    
15. class FFT
    
16. {
    
17. private:
    
18. void changeOrder(double* xr,double* xi,int n);
    
19. public:
    
20. FFT(void);
    
21. void FFT_1D(double* ctxr,double* ctxi,double* cfxr,double* cfxi,int len);
    
22. void rFFT_1D(double* ctxr,double* cfxr,double* cfxi,int len);
    
23. void IFFT_1D(double* cfxr,double* cfxi,double* ctxr,double* ctxi,int len);
    
24. void rIFFT_1D(double* cfxr,double* cfxi,double* ctxr,int len);
    
25. ~FFT(void);
    
26. };
    
27. #endif
    
28. 2.FFT.cpp
    
29. #include ".fft.h"
    
30. FFT::FFT()
    
31. {
    
32. }
    
33. ///////////////////////////
    
34. //倒序实现
    
35. //xr实部，xi虚部，n为2的幂
    
36. ///////////////////////////
    
37. void FFT::changeOrder(double *xr,double *xi,int n)
    
38. {
    
39. double temp;
    
40. int k;
    
41. int lh=n/2;
    
42. int j=lh;
    
43. int n1=n-2;
    
44. for(int i=1;i<=n1;i++)
    
45. {
    
46.   if(i<j)
    
47.   {
    
48.    temp=xr[i];
    
49.    xr[i]=xr[j];
    
50.    xr[j]=temp;
    
51.    temp=xi[i];
    
52.    xi[i]=xi[j];
    
53.    xi[j]=temp;
    
54.   }
    
55.   k=lh;
    
56.   while(j>=k)
    
57.   {
    
58.    j=j-k;
    
59.    k=(int)(k/2+0.5);
    
60.   }
    
61.   j=j+k;
    
62. }
    
63. }
    
64. //////////////////////////
    
65. //复数FFT
    
66. //ctxr和ctxi的长度为len
    
67. //cfxr和cfxi的长度为2的幂
    
68. //////////////////////////
    
69. void FFT::FFT_1D(double *ctxr,double *ctxi,double *cfxr,double *cfxi,int len)
    
70. {
    
71. int m=ceil(log((double)len)/log(2.0));
    
72. int l,b,j,p,k;
    
73. double rkb,ikb;
    
74. int n=1<<m;
    
75. double* rcos=new double[n/2];
    
76. double* isin=new double[n/2];
    
77. for(l=0;l<n/2;l++)           
    
78. {
    
79.   rcos[l]=cos(l*pi*2/n);
    
80.   isin[l]=sin(l*pi*2/n);
    
81. }                           
    
82. memcpy(cfxr,ctxr,sizeof(double)*len);
    
83. memcpy(cfxi,ctxi,sizeof(double)*len);
    
84. for(l=len;l<n;l++)
    
85. {
    
86.   cfxr[l]=0;
    
87.   cfxi[l]=0;
    
88. }
    
89. changeOrder(cfxr,cfxi,n);  //倒序
    
90. for(l=1;l<=m;l++)
    
91. {
    
92.   b=(int)(pow(2,l-1)+0.5);
    
93.   for(j=0;j<b;j++)
    
94.   {
    
95.    p=j*(int)(pow(2,m-l)+0.5);
    
96.    for(k=j;k<n;k+=(int)(pow(2,l)+0.5))
    
97.    {
    
98.     rkb=cfxr[k+b]*rcos[p]+cfxi[k+b]*isin[p];
    
99.     ikb=cfxi[k+b]*rcos[p]-cfxr[k+b]*isin[p];
    
100.     cfxr[k+b]=cfxr[k]-rkb;
     
101.     cfxi[k+b]=cfxi[k]-ikb;
     
102.     cfxr[k]=cfxr[k]+rkb;
     
103.     cfxi[k]=cfxi[k]+ikb;
     
104.    }
     
105.   }
     
106. }
     
107. delete []rcos;
     
108. delete []isin;
     
109. }
     
110. /////////////////////////////
     
111. //实数FFT
     
112. //ctxr的长度为len
     
113. //cfxr和cfxi的长度为2的幂
     
114. ////////////////////////////
     
115. void FFT::rFFT_1D(double *ctxr,double *cfxr,double *cfxi,int len)
     
116. {
     
117. int m=ceil(log((double)len)/log(2.0));
     
118. int n=1<<m;
     
119. double* rcos=new double[n/2];
     
120. double* isin=new double[n/2];
     
121. for(int l=0;l<n/2;l++)         
     
122. {
     
123.   rcos[l]=cos(l*pi*2/n);
     
124.   isin[l]=sin(l*pi*2/n);
     
125. }   
     
126. double* txr=new double[(len+1)/2];
     
127. double* txi=new double[(len+1)/2];
     
128. for(int i=0;i<len/2;i++)
     
129. {
     
130.   txr[i]=ctxr[i*2];
     
131.   txi[i]=ctxr[i*2+1];
     
132. }
     
133. if(len%2==1)
     
134. {
     
135.   txr[(len-1)/2]=ctxr[len-1];
     
136.   txi[(len-1)/2]=0;
     
137. }
     
138. FFT_1D(txr,txi,cfxr,cfxi,(len+1)/2);
     
139. double* x1r=new double[n/2];
     
140. double* x1i=new double[n/2];
     
141. double* x2r=new double[n/2];
     
142. double* x2i=new double[n/2];
     
143. x1r[0]=cfxr[0];
     
144. x1i[0]=0;
     
145. x2r[0]=cfxi[0];
     
146. x2i[0]=0;
     
147. for(int k=1;k<n/2;k++)
     
148. {
     
149.   x1r[k]=(cfxr[k]+cfxr[n/2-k])/2.0;
     
150.   x1i[k]=(cfxi[k]-cfxi[n/2-k])/2.0;
     
151.   x2r[k]=(cfxi[k]+cfxi[n/2-k])/2.0;
     
152.   x2i[k]=(-cfxr[k]+cfxr[n/2-k])/2.0;
     
153. }
     
154. double rkb,ikb;
     
155. for(i=0;i<n/2;i++)
     
156. {
     
157.   rkb=x2r[i]*rcos[i]+x2i[i]*isin[i];
     
158.   ikb=x2i[i]*rcos[i]-x2r[i]*isin[i];
     
159.   cfxr[i+n/2]=x1r[i]-rkb;
     
160.   cfxi[i+n/2]=x1i[i]-ikb;
     
161.   cfxr[i]=x1r[i]+rkb;
     
162.   cfxi[i]=x1i[i]+ikb;
     
163. }
     
164. delete []txr;
     
165. delete []txi;
     
166. delete []x1r;
     
167. delete []x1i;
     
168. delete []x2r;
     
169. delete []x2i;
     
170. delete []rcos;
     
171. delete []isin;
     
172. }
     
173. ///////////////////////////////
     
174. //复数IFFT
     
175. //cfxr和cfxi的长度为n(2的幂)
     
176. //ctxr和ctxi的长度为len
     
177. ///////////////////////////////
     
178. void FFT::IFFT_1D(double *cfxr,double *cfxi,double *ctxr,double *ctxi,int len)
     
179. {
     
180. int m=ceil(log((double)len)/log(2.0));
     
181. int n=1<<m;
     
182. double* txr=new double[n];
     
183. double* txi=new double[n];
     
184. for(int i=0;i<n;i++)
     
185.   cfxi[i]=-cfxi[i];
     
186. FFT_1D(cfxr,cfxi,txr,txi,n);
     
187. for(i=0;i<len;i++)
     
188. {
     
189.   ctxr[i]=txr[i]/n;
     
190.   ctxi[i]=-txi[i]/n;
     
191. }
     
192. delete []txr;
     
193. delete []txi;
     
194. }
     
195. //////////////////////////////
     
196. //实数IFFT
     
197. //cfxr和cfxi的长度为n(2的幂)
     
198. //ctxr的长度为len
     
199. //////////////////////////////
     
200. void FFT::rIFFT_1D(double *cfxr,double *cfxi,double *ctxr,int len)
     
201. {
     
202. int m=ceil(log((double)len)/log(2.0));
     
203. int n=1<<m;
     
204. double* txr=new double[n];
     
205. double* txi=new double[n];
     
206. for(int i=0;i<n;i++)
     
207.   cfxi[i]=-cfxi[i];
     
208. FFT_1D(cfxr,cfxi,txr,txi,n);
     
209. for(int i=0;i<len;i++)
     
210. {
     
211.   ctxr[i]=txr[i]/n;
     
212. }
     
213. delete []txr;
     
214. delete []txi;
     
215. }
     
216. FFT::~FFT(void)
     
217. {
     
218. }
     
219. 3.test.cpp
     
220. #include <iostream>
     
221. #include "FFT.h"
     
222. using namespace std;
     
223. void main()
     
224. {
     
225. double xr[10]={1,2,3,4,5,6,7,8,9,10};   //实部
     
226. double xi[10]={0,0,0,0,0,0,0,0,0,0};    //虚部
     
227. double *cxr,*cxi;
     
228. cxr=new double[16];
     
229. cxi=new double[16];
     
230. FFT f;
     
231. f.rFFT_1D(xr,cxr,cxi,10);
     
232. for(int i=0;i<16;i++)
     
233. {
     
234.   cout<<cxr[i]<<"+j"<<cxi[i];
     
235.   cout<<endl;
     
236. }
     
237. cout<<endl;
     
238. double *fxr,*fxi;
     
239. fxr=new double[10];
     
240. fxi=new double[10];
     
241. f.rIFFT_1D(cxr,cxi,fxr,10);
     
242. for(i=0;i<10;i++)
     
243. {
     
244.   cout<<fxr[i];
     
245.   cout<<endl;
     
246. }
     
247. delete []fxr;
     
248. delete []fxi;
     
249. delete []cxr;
     
250. delete []cxi;
     
251. 
     
252. }

复制代码

lebronze

多谢各位的热心回复，很感动。
不过2的整数幂的fft我也会做，现在的问题是如果信号长度是任意的，可能不是2的整数幂，那要
怎么做fft？？
matlab的fft函数就可以实现这个功能，但它内部是封装好的看不到代码。我查网上说是通过多种
混合基（不同长度的基混合）来做的fft。
希望有经验的人能指点迷津

lebronze

eastar 发表于 2016-7-6 08:18
参考fft matlab代码改成c行吗

matlab的fft是封装好的，看不到啊

lebronze

impulse 发表于 2016-7-6 12:36
matlab fft用的是FFTW

好的，那我还是下个FFTW看看他处理非2整数幂fft的思路吧

impulse

lebronze 发表于 2016-7-6 15:19
好的，那我还是下个FFTW看看他处理非2整数幂fft的思路吧

FFTW好像没有嵌入式版本

eastar

lebronze 发表于 2016-7-6 15:18
matlab的fft是封装好的，看不到啊

我看论坛里  之前有人问fft源代码的

Generation

你可以参考一下张明的《任意长度FFT算法的C++实现》

lebronze

Generation 发表于 2016-7-7 08:46
你可以参考一下张明的《任意长度FFT算法的C++实现》

多谢回复，我这两天看了一下张明的代码。
他相当于使用多组基的fft来做，但是问题是它的代码是C++的，我本身就看不太懂，在自己用C写真的很困难。。。
我在想可不可以直接用DFT来做