声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 5024|回复: 52

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

[复制链接]
发表于 2016-7-5 21:56 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

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

本帖被以下淘专辑推荐:

回复
分享到:

使用道具 举报

 楼主| 发表于 2016-7-13 10:21 | 显示全部楼层
多谢各位回复,目前的解决方法是:对原始信号末尾补零到2的整数幂,在进行FFT。
时域补零,相当于频域插值,不会改变频谱的包络

点评

明白就好。  发表于 2016-7-13 21:06
发表于 2016-7-5 22:34 | 显示全部楼层
本帖最后由 impulse 于 2016-7-5 22:38 编辑

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

多谢回复,不过我这边是要搞嵌入式开发,FFT,IFFT只能自己写,不能调用库
发表于 2016-7-6 07:56 | 显示全部楼层
从编程角度看,FFT,IFFT 的精髓就是必须是2的整数幂,加零等方法都是凑数的近似方法。
发表于 2016-7-6 08:18 | 显示全部楼层
参考fft matlab代码改成c行吗

点评

有现成的C语言FFT代码。  发表于 2016-7-6 11:11
发表于 2016-7-6 12:36 | 显示全部楼层
eastar 发表于 2016-7-6 08:18
参考fft matlab代码改成c行吗

matlab fft用的是FFTW
发表于 2016-7-6 12:56 | 显示全部楼层
发表于 2016-7-6 13:11 | 显示全部楼层
这是我找到的C++的代码,不知道有没有用  你看一下
  1. 时域抽取法FFT的C++实现。
  2.      下面是算法的流程图

  3.      倒序的流程图

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

  249. }
复制代码

评分

1

查看全部评分

 楼主| 发表于 2016-7-6 15:17 | 显示全部楼层
多谢各位的热心回复,很感动。
不过2的整数幂的fft我也会做,现在的问题是如果信号长度是任意的,可能不是2的整数幂,那要
怎么做fft??
matlab的fft函数就可以实现这个功能,但它内部是封装好的看不到代码。我查网上说是通过多种
混合基(不同长度的基混合)来做的fft。
希望有经验的人能指点迷津

点评

如果你想完成任意长度的数据是很难实现的,Matlab中关于这段代码从来不显示。我常用的方法就是加零数据延长到2的整数次幂。你也可以试一下,效果还可以!!!  详情 回复 发表于 2016-7-12 11:18
如果要计较计算速度的话,补零是最佳的途径。  发表于 2016-7-9 09:09
 楼主| 发表于 2016-7-6 15:18 | 显示全部楼层
eastar 发表于 2016-7-6 08:18
参考fft matlab代码改成c行吗

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

点评

我看论坛里 之前有人问fft源代码的  详情 回复 发表于 2016-7-7 08:38
 楼主| 发表于 2016-7-6 15:19 | 显示全部楼层
impulse 发表于 2016-7-6 12:36
matlab fft用的是FFTW

好的,那我还是下个FFTW看看他处理非2整数幂fft的思路吧
发表于 2016-7-6 20:03 | 显示全部楼层
lebronze 发表于 2016-7-6 15:19
好的,那我还是下个FFTW看看他处理非2整数幂fft的思路吧

FFTW好像没有嵌入式版本
发表于 2016-7-7 08:38 | 显示全部楼层
lebronze 发表于 2016-7-6 15:18
matlab的fft是封装好的,看不到啊

我看论坛里  之前有人问fft源代码的
发表于 2016-7-7 08:46 | 显示全部楼层
你可以参考一下张明的《任意长度FFT算法的C++实现》
 楼主| 发表于 2016-7-8 22:36 | 显示全部楼层
Generation 发表于 2016-7-7 08:46
你可以参考一下张明的《任意长度FFT算法的C++实现》

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

点评

我有个同学就把matlab程序改为c++的了 他说学习起来还不算太难  详情 回复 发表于 2016-7-11 11:24
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-15 12:55 , Processed in 0.118179 second(s), 37 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表