fft

最新推荐文章于 2023-02-21 17:55:57 发布
最新推荐文章于 2023-02-21 17:55:57 发布
阅读量1k 收藏 1
点赞数
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/q123456789098/article/details/50605890
版权 
#include <stdio.h> 
#include <stdlib.h> 
#include <math.h> 
//const float PI = 3.1416; 


inline void swap (float &a, float &b) 
{ 
	float t; 
	t = a; 
	a = b; 
	b = t; 
} 
void bitrp (float xreal [], float ximag [], int n) 
{ 
	//位反转置换Bit-reversal Permutation 
	int i, j, a, b, p; 
	for (i = 1, p = 0; i < n; i *= 2) 
	{ 
		p ++;
	} 
	for (i = 0; i < n; i ++) 
	{ 
		a = i; 
		b = 0; 
		for (j = 0; j < p; j ++) 
		{ 
			b = (b << 1) + (a & 1); // b = b * 2 + a % 2; 
			a >>= 1; // a = a / 2; 
		} 
		if ( b > i) 
		{ 
			swap (xreal [i], xreal [b]); 
			swap (ximag [i], ximag [b]); 
		} 
	} 
}

void FFT(float xreal[], float ximag[], int n) 
{ 
	//快速傅立叶变换，将复数x 变换后仍保存在x 中，xreal, ximag 分别是x 的实部和虚部
	float wreal[128];
	float wimag[128];
	float treal,timag,ureal,uimag,arg;
	int m, k, j, t, index1, index2; 
	bitrp (xreal, ximag, n); 
	// 计算1 的前n / 2 个n 次方根的共轭复数W'j = wreal [j] + i * wimag [j] , j = 0, 1, , n / 2 - 1 
	
	arg = - 2 * PI / n; 
	treal = cos (arg); 
	timag = sin (arg); 
	wreal [0] = 1.0; 
	wimag [0] = 0.0; 
	for (j = 1; j < n / 2; j ++) 
	{ 
		wreal [j] = wreal [j - 1] * treal - wimag [j - 1] * timag; 
		wimag [j] = wreal [j - 1] * timag + wimag [j - 1] * treal; 
	} 
	for (m = 2; m <= n; m *= 2)  
	{ 
		for (k = 0; k < n; k += m)
		{ 
			for (j = 0; j < m / 2; j ++) 
			{ 
				index1 = k + j;
				index2 = index1 + m / 2; 
				t = n * j / m; 
				// 旋转因子w 的实部在wreal [] 中的下标为t 
				
				treal = wreal [t] * xreal [index2] - wimag [t] * ximag [index2]; 
				timag = wreal [t] * ximag [index2] + wimag [t] * xreal [index2]; 
				ureal = xreal [index1]; 
				uimag = ximag [index1]; 
				xreal [index1] = ureal + treal; 
				ximag [index1] = uimag + timag; 
				xreal [index2] = ureal - treal; 
				ximag [index2] = uimag - timag; 
			} 
		} 
	} 
} 

// 
// void IFFT (float xreal [], float ximag [], int n) 
// { 
// 	// 快速傅立叶逆变换
// 	float wreal [128], wimag [128], treal, timag, ureal, uimag, arg; 
// 	int m, k, j, t, index1, index2; 
// 	bitrp (xreal, ximag, n); 
// 	// 计算1 的前n / 2 个n 次方根Wj = wreal [j] + i * wimag [j] , j = 0, 1, , n / 2 - 1 
// 	
// 	arg = 2 * PI / n; 
// 	treal = cos (arg); 
// 	timag = sin (arg); 
// 	wreal [0] = 1.0; 
// 	wimag [0] = 0.0; 
// 	for (j = 1; j < n / 2; j ++) 
// 	{ 
// 		wreal [j] = wreal [j - 1] * treal - wimag [j - 1] * timag; 
// 		wimag [j] = wreal [j - 1] * timag + wimag [j - 1] * treal; 
// 	} 
// 	for (m = 2; m <= n; m *= 2) 
// 	{ 
// 		for (k = 0; k < n; k += m)  
// 		{ 
// 			for (j = 0; j < m / 2; j ++)
// 			{ 
// 				index1 = k + j; 
// 				index2 = index1 + m / 2; 
// 				t = n * j / m; 
// 				// 旋转因子w 的实部在wreal [] 中的下标为t 
// 				treal = wreal [t] * xreal [index2] - wimag [t] * ximag [index2];
// 				timag = wreal [t] * ximag [index2] + wimag [t] * xreal [index2]; 
// 				ureal = xreal [index1]; 
// 				uimag = ximag [index1];
// 				xreal [index1] = ureal + treal;
// 				ximag [index1] = uimag + timag; 
// 				xreal [index2] = ureal - treal; 
// 				ximag [index2] = uimag - timag; 
// 			} 
// 		} 
// 	} 
// 	for (j=0; j < n; j ++) 
// 	{ 
// 		xreal [j] /= n; 
// 		ximag [j] /= n; 
// 	} 
// } 


    #include "math.h"  
    #include "fft.h"  
    //精度0.0001弧度  
      
    void conjugate_complex(int n,complex in[],complex out[])  
    {  
      int i = 0;  
      for(i=0;i<n;i++)  
      {  
        out[i].imag = -in[i].imag;  
        out[i].real = in[i].real;  
      }   
    }  
      
    void c_abs(complex f[],float out[],int n)  
    {  
      int i = 0;  
      float t;  
      for(i=0;i<n;i++)  
      {  
        t = f[i].real * f[i].real + f[i].imag * f[i].imag;  
        out[i] = sqrt(t);  
      }   
    }  
      
      
    void c_plus(complex a,complex b,complex *c)  
    {  
      c->real = a.real + b.real;  
      c->imag = a.imag + b.imag;  
    }  
      
    void c_sub(complex a,complex b,complex *c)  
    {  
      c->real = a.real - b.real;  
      c->imag = a.imag - b.imag;   
    }  
      
    void c_mul(complex a,complex b,complex *c)  
    {  
      c->real = a.real * b.real - a.imag * b.imag;  
      c->imag = a.real * b.imag + a.imag * b.real;     
    }  
      
    void c_div(complex a,complex b,complex *c)  
    {  
      c->real = (a.real * b.real + a.imag * b.imag)/(b.real * b.real +b.imag * b.imag);  
      c->imag = (a.imag * b.real - a.real * b.imag)/(b.real * b.real +b.imag * b.imag);  
    }  
      
    #define SWAP(a,b)  tempr=(a);(a)=(b);(b)=tempr  
      
    void Wn_i(int n,int i,complex *Wn,char flag)  
    {  
      Wn->real = cos(2*PI*i/n);  
      if(flag == 1)  
      Wn->imag = -sin(2*PI*i/n);  
      else if(flag == 0)  
      Wn->imag = -sin(2*PI*i/n);  
    }  
      
    //傅里叶变化  
    void fft(int N,complex f[])  
    {  
      complex t,wn;//中间变量  
      int i,j,k,m,n,l,r,M;  
      int la,lb,lc;  
      /*----计算分解的级数M=log2(N)----*/  
      for(i=N,M=1;(i=i/2)!=1;M++);   
      /*----按照倒位序重新排列原信号----*/  
      for(i=1,j=N/2;i<=N-2;i++)  
      {  
        if(i<j)  
        {  
          t=f[j];  
          f[j]=f[i];  
          f[i]=t;  
        }  
        k=N/2;  
        while(k<=j)  
        {  
          j=j-k;  
          k=k/2;  
        }  
        j=j+k;  
      }  
      
      /*----FFT算法----*/  
      for(m=1;m<=M;m++)  
      {  
        la=pow(2,m); //la=2^m代表第m级每个分组所含节点数       
        lb=la/2;    //lb代表第m级每个分组所含碟形单元数  
                     //同时它也表示每个碟形单元上下节点之间的距离  
        /*----碟形运算----*/  
        for(l=1;l<=lb;l++)  
        {  
          r=(l-1)*pow(2,M-m);     
          for(n=l-1;n<N-1;n=n+la) //遍历每个分组，分组总数为N/la  
          {  
            lc=n+lb;  //n,lc分别代表一个碟形单元的上、下节点编号       
            Wn_i(N,r,&wn,1);//wn=Wnr  
            c_mul(f[lc],wn,&t);//t = f[lc] * wn复数运算  
            c_sub(f[n],t,&(f[lc]));//f[lc] = f[n] - f[lc] * Wnr  
            c_plus(f[n],t,&(f[n]));//f[n] = f[n] + f[lc] * Wnr  
          }  
        }  
      }  
    }  
      
    //傅里叶逆变换  
    void ifft(int N,complex f[])  
    {  
      int i=0;  
      conjugate_complex(N,f,f);  
      fft(N,f);  
      conjugate_complex(N,f,f);  
      for(i=0;i<N;i++)  
      {  
        f[i].imag = (f[i].imag)/N;  
        f[i].real = (f[i].real)/N;  
      }  
    }

matlab仿真

%例1：x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)。采样频率fs=100Hz，分别绘制N=128、1024点幅频图。

clf;
fs=100;N=128;   %采样频率和数据点数
n=0:N-1;t=n/fs;   %时间序列
x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号
y=fft(x,N);    %对信号进行快速Fourier变换
mag=abs(y);     %求得Fourier变换后的振幅
f=n*fs/N;    %频率序列
subplot(2,2,1),plot(f,mag);   %绘出随频率变化的振幅
xlabel('频率/Hz');
ylabel('振幅');title('N=128');grid on;
subplot(2,2,2),plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅
xlabel('频率/Hz');
ylabel('振幅');title('N=128');grid on;
%对信号采样数据为1024点的处理
fs=100;N=1024;n=0:N-1;t=n/fs;
x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号
y=fft(x,N);   %对信号进行快速Fourier变换
mag=abs(y);   %求取Fourier变换的振幅
f=n*fs/N;
subplot(2,2,3),plot(f,mag); %绘出随频率变化的振幅
xlabel('频率/Hz');
ylabel('振幅');title('N=1024');grid on;
subplot(2,2,4)
plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅
xlabel('频率/Hz');
ylabel('振幅');title('N=1024');grid on;


DFT

#include <stdio.h>  
#include <math.h> 
#include "fft.h" 
#define N   SAMPLEPOINTS
              //序列长度   
 
    
typedef float ElementType; //原始数据序列的数据类型,可以在这里设置   
typedef struct              //复数结构体,用于实现傅里叶运算   
{  
    ElementType real,imag;  
}complex;  
complex dataResult[N];      //傅里叶运算的频域结果序列的值(复数)  
ElementType dataSource[N];  //输入的原始数据序列   
ElementType dataFinualResult[N]; //最终频域序列结果(dataResult[N]的模,实数)  
void FFT_Calculate_OneNode(int k)//计算频域上一个点的DFT值   
{  
    int n = 0;  
    complex ResultThisNode;  
    complex part[N];  
    ResultThisNode.real = 0;  
    ResultThisNode.imag = 0;  
    for(n=0; n<N; n++)  
    {  
        part[n].real = cos(2*PI/N*k*n)*dataSource[n];//运用欧拉公式把复数拆分成实部和虚部   
        part[n].imag = sin(2*PI/N*k*n)*dataSource[n];  
		
        ResultThisNode.real += part[n].real;  
        ResultThisNode.imag += part[n].imag;  
    }  
    dataResult[k].real = ResultThisNode.real;  
    dataResult[k].imag = ResultThisNode.imag;  
}  
void FFT_Calculate()//对输入序列全部计算DFT  
{  
    int i = 0;  
    for(i=0; i<N; i++)  
    {  
        FFT_Calculate_OneNode(i);  
        dataFinualResult[i] = sqrt(dataResult[i].real * dataResult[i].real + dataResult[i].imag * dataResult[i].imag);  
    }  
}  
int testdftmain(float * input,float * output)  
{  
    int i = 0;  
    for(i=0; i<N; i++)//制造输入序列   
    {  
        dataSource[i] = input[i];  
        
    }  
    FFT_Calculate();//进行DFT计算   
    
    for(i=0; i<N; i++)  
	{
		output[i]=dataFinualResult[i];
	}  
    return 0;  
}



q123456789098
关注
  • 0
    点赞
  • 踩
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
十分简明易懂的FFT（快速傅里叶变换）
路人黑的纸巾
08-07 59万+
FFT有什么用 快速傅里叶变换 (fast Fourier transform),即利用计算机计算离散傅里叶变换（DFT)的高效、快速计算方法的统称，简称FFT。快速傅里叶变换是1965年由J.W.库利和T.W.图基提出的。采用这种算法能使计算机计算离散傅里叶变换所需要的乘法次数大为减少，特别是被变换的抽样点数N越多，FFT算法计算量的节省就越显著。 FFT（Fast Fouri...
MATLAB仿真文件，对IQ数据进行FFT
03-11
在本项目中，"MATLAB仿真文件，对IQ数据进行FFT" 提供了一组用于处理IQ（In-phase and Quadrature）数据的代码和模型，通过快速傅里叶变换(FFT)分析信号的频域特性。 首先，我们来理解一下IQ数据。在通信系统中，IQ...
JAVA实现FFT算法
01-03
java实现FFT算法，关于快速傅里叶变换（FFT）和傅里叶变换的理论知识这里我就不提了，本文主要讲解FFT实现： 之前想找一个FFT代码，在网上找了很多都是有问题的，下面我完善了一个供大家学习交流；
FFT的详细解释，相信你看了就明白了
teresa_zp的专栏
12-28 11万+
FFT是离散傅立叶变换的快速算法，可以将一个信号变换到频域。有些信号在时域上是很难看出什么特征的，但是如果变换到频域之后，就很容易看出特征了。这就是很多信号分析采用FFT变换的原因。另外，FFT可以将一个信号的频谱提取出来，这在频谱分析方面也是经常用的。 虽然很多人都知道F.........
FFT（快速傅里叶变换）
weixin_42638946的博客
10-17 1万+
FFT 1. FFT原理 原理 FFT（Fast Fourier Transformation），中文名快速傅里叶变换，用来加速多项式乘法。 朴素高精度乘法时间复杂度是O(n2)O(n ^ 2)O(n2)，n是数据位数；但FFT能在O(n×log(n))O(n \times log(n))O(n×log(n))的时间内解决。 假设给定一个多项式： A(x)=a0+a1x+a2x2+...an−1xn−1 A(x) = a_0 + a_1 x + a_2 x^2 + ... a_{n-1}
FFT简介小结
知识搬运工的博客
04-18 9040
FFT是离散傅立叶变换的快速算法，可以将一个信号变换到频域。有些信号在时域上是很难看出什么特征的，但是如果变换到频域之后，就很容易看出特征了。这就是很多信号分析采用FFT变换的原因。另外，FFT可以将一个信号的频谱提取出来，这在频谱分析方面也是经常用的。虽然很多人都知道FFT是什么，可以用来做什么，怎么去做，但是却不知道FFT之后的结果是什意思、如何决定要使用多少点来做FFT。 ...
快速傅里叶变换——FFT
upupyon996deqing的博客
09-13 2万+
2，用一个公式来表示一下，这个公式的意义就是用来分解信号的，因为之前说FFT使用一个个正弦波来表示时域信号的，所以我们要确定正弦信号和时域信号的相关性，并且这个正弦波也要匹配块长度。上面的这个图就是加上传函数之后出现的FFT频谱图，但是可以看到激波的Bin比较宽，出现这个线性的原因是因为和快长度有关，如果提高块长度，那么分辨率就会变高，那么就会出现一个比较窄的Bin。其中，块长度也就是样本的数量跟FFT采集的周期有直接的关系，块长度对频率分辨率有直接的影响，假设采样的频率是48kHZ,块长度是512。
FFT.rar_fft
09-21
快速傅立叶变换（FFT）是一种高效的计算离散傅立叶变换（DFT）和其逆变换的方法，广泛应用于数字信号处理、图像处理、通信工程等多个领域。在本压缩包"FFT.rar_fft"中，包含了在TI公司的TMS320C6713数字信号处理器上...
FFT_fft_matlab_
09-30
标题中的"FFT_fft_matlab_"表明我们关注的主题是快速傅里叶变换（FFT）在MATLAB环境中的应用以及如何用C语言实现相同的功能。FFT是一种高效计算离散傅里叶变换（DFT）的算法，广泛应用于信号处理、图像处理、数字...
STM32的FFT官方库.zip
04-29
例如，可能有`fft_init()`用于设置FFT配置，`fft_execute()`用于执行变换，`fft_result()`用于获取变换结果等。 5. **性能优化**：STM32的FFT库可能已经针对特定的处理器进行了优化，但开发者还可以根据具体应用...
Test_FFT_fft_
10-01
标题中的"Test_FFT_fft_"暗示我们关注的主题是快速傅里叶变换（FFT），这是一个在数字信号处理领域中至关重要的算法。FFT是一种高效计算离散傅里叶变换（DFT）的方法，它大大减少了计算复杂度，使得大规模数据的频谱...
FFT算法讲解——麻麻我终于会FFT了！
热门推荐
谁能代表肯德基
03-12 15万+
FFT——快速傅里叶变换 这块不写东西空荡荡的，我决定还是把FFT的定义给贴上吧 FFT（Fast Fourier Transformation）是离散傅氏变换（DFT）的快速算法。即为快速傅氏变换。它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。 这三段话其实一点用也没有 FFT是干什么的 FFT在算法竞赛中就有一个用途：加速多项式乘法 ...
FFT(傅里叶快速变换，详细讲解+推导) 每日一遍，算法再见！
TheWayForDream
02-22 2万+
FFT详细推导FFT(傅里叶快速变换)一.前置知识1.复数和单位根2.单位根的三个引理3.多项式二.FFT(快速傅里叶变换推导)三.IFFT四.FFT求解多项式乘积模板代码1.递归版2.非递归版(这个更快，省去了递归时间)五.视频资源 FFT(傅里叶快速变换) FFT在实际工程中有着非常的广泛，尤其是在信号领域，在ACM算法竞赛领域主要可以用来快速计算多项式的乘积 一.前置知识 1.复数和单位根 有人会觉得FFT怎么会用到复数的知识，想必是有点古怪，后面在推导过程中会介绍到它的用处。 2.单位根的三个引
快速傅里叶变换 （FFT）基础
weixin_46118768的博客
07-15 1万+
FFT
一文让你秒懂FFT（快速傅里叶变换）
S2191300319的博客
02-21 8682
一文带你领略傅里叶变换的魅力
FFT频谱分析
蜗牛在听雨的博客
05-24 4万+
FFT是离散傅立叶变换的快速算法，可以将一个信号变换到频域。有些信号在时域上是很难看出什么特征的，但是如果变换到频域之后，就很容易看出特征了。这就是很多信号分析采用FFT变换的原因。另外，FFT可以将一个信号的频谱提取出来，这在频谱分析方面也是经常用的。 虽然很多人都知道FFT是什么，可以用来做什么，怎么去做，但是却不知道FFT之后的结果是什意思、如何决定要使用多少点来做FFT。 现在就根据实际经...
计算机模拟考试.doc
最新发布
07-18
计算机
快速傅里叶变换(FFT)：DIT-FFT、DIF-FFT与应用
"本资源主要介绍了快速傅里叶变换（FFT）的概念，包括DIT-FFT、DIF-FFT、IFFT、Chirp-FFT和线性卷积的FFT算法，并探讨了FFT在信号处理中的应用，如频谱分析、滤波器实现和实时信号处理。" 快速傅里叶变换（FFT）是...

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值