Java快速傅里叶变换FFT的程序实现（时间抽取的基-2FFT倒位计算蝶形运算）

Posted 2022-05-18 阿汪先生

Java——快速傅里叶变换（FFT）的程序实现

---

好久没来更新了，阿汪大三了。
这学期阿汪要学习两门课《数字信号处理》和《Java程序设计》，刚好前几天老师告诉我们不久后会有个实验，要求我们编写一个程序实现快速傅里叶变换（FFT）。
所以，阿汪用Java写了一下。
//废话和原理就不多说了，直接上程序！！！

---

//阿汪先生的博客.ws
import java.util.Scanner;
import java.math.*;

//对数函数
 class Log
		   static public double log(double value, double base) 
		return Math.log(value) / Math.log(base);
		   
 
 
//乘方运算
class Pow	
	 public static int powNM(int n,int m)
		    if(m == 1)
		        return n;
		    else
		        return n*powNM(n,m-1);
		    //Math.pow(m,n);
		


public class DFT_FFT 	
	public static void main (String []args)
		//初始化变量
		final double PI=3.141;
		int N1=0;	//序列的长度
		int M=0;   //序列的级数，数组的长度为2^M
		int arraylength; //数组长度
		int i=0;	//输入数组的循环次数
		int t=0;  //输出数组的循环次数						

		
///	
		
		//输入合法化
		Scanner sc=new Scanner (System.in);
		while(N1<2)
			System.out.println("请确定你输入序列的长度！");
			N1=sc.nextInt();
				
///	
		
		//判断长度并初始化数组
		//确定N=2^M的级数M
		M=(Log.log(N1, 2)==(int) Log.log(N1, 2))?(int) Log.log(N1, 2) : ((int)Log.log(N1, 2))+1;
		//确定数字长度
		arraylength=(int)Pow.powNM(2, M);							
		System.out.printf("序列长度N1为%d时，数组长度为%d，级数M为%d !\\n",N1,arraylength,M);
	    float []Array=new  float [ arraylength ];   //实部
		float []IArray=new  float [ arraylength ];  //虚部
///	
       //输入序列
		float jc;   //暂存转换后的数据
		String xy;
	    System.out.printf("请依次输入序列的%d个值：\\n",N1);
	    System.out.println("");
		while(i<N1)
			System.out.printf("x[%d]实部：\\n",i);
			xy=sc.next();
			if((xy.charAt(0)<48|xy.charAt(0)>57))   //输入不是数字
				System.out.println("阿汪先生(a_wang_xian_sheng)的博客！");
				break;			//跳出循环输入
			
			if((xy.charAt(0)>=48|xy.charAt(0)<=57)) //输入数字
				jc = Float.parseFloat(xy);			//String转化成float类型
				Array[i]=jc;						//保存String转存下来的数字
			
			System.out.printf("x[%d]虚部：\\n",i);   //输入虚部
			IArray[i]=sc.nextFloat();				//录入数字
				i++;
			
			
		
		  System.out.println("");
		  System.out.println("输入结果如下：");
			while(t<arraylength)
				 System.out.printf("x[%d]=",t);
				 System.out.print(Array[t]);
				 System.out.printf("+j(%f)\\n",IArray[t]);
				 t++;
			
//	
			int w=0;   //倒序处理时的外循环次数
			int j=0;   //倒序处理时的内循环次数
			float ws=0;  //倒序处理时的数组间值传递变量
			int p=0;   //倒序处理时下一倒位序的确定变量
			int w1=0;   //倒序处理时的外循环次数
			int j1=0;   //倒序处理时的内循环次数
			float ws1=0;  //倒序处理时的数组间值传递变量
			int p1=0;   //倒序处理时下一倒位序的确定变量
		//倒序处理(实部)
	  for(w=1,j=arraylength/2;w<arraylength-1;w++)            
	   
	     if(w<j)                    //如果i<j,即进行变址    
	     
	       ws=Array[j];                  
	       Array[j]=Array[w];      
	       Array[w]=ws;    
	     
	       p=arraylength/2;                    //求j的下一个倒位序 
	     while(p<=j)               //如果k<=j,表示j的最高位为1      
	               
	       j=j-p;                 //把最高位变成0 
	       p=p/2;                 //k/2，比较次高位，依次类推，逐个比较，直到某个位为0      
	   
	    j=j+p;                   //把0改为1  
	  
		
		//倒序处理(虚部)
	  for(w1=1,j1=arraylength/2;w1<arraylength-1;w1++)            
	   
	     if(w1<j1)                    //如果i<j,即进行变址    
	     
	       ws1=IArray[j1];                  
	       IArray[j1]=IArray[w1];      
	       IArray[w1]=ws1;    
	     
	       p1=arraylength/2;                    //求j的下一个倒位序 
	     while(p1<=j1)               //如果k<=j,表示j的最高位为1      
	               
	       j1=j1-p1;                 //把最高位变成0 
	       p1=p1/2;                 //k/2，比较次高位，依次类推，逐个比较，直到某个位为0      
	   
	    j1=j1+p1;                   //把0改为1  
	  
//
		
	  //蝶形运算变量定义并初始化
		int L=1;//蝶形运算级数，用于循环
		int N=2;//蝶形运算数据量，用于循环 ，WN(k)*X(k+N/2)中的N
		int distance=1;//蝶形运算两节点间的距离，用于循环（distance=N/2） 
		int group=arraylength/2;  //蝶形运算的组数，长度为数组长度的一半
		double tempr1,tempr2,tempi1,tempi2;//临时变量
		int k=0;  //WN(k)*X(k+N/2)中的k	
	 	
		
	  //蝶形运算
		for(;L<=M;L++)//第L级的蝶形运算
			for(i=0;i<group;i++)//第group组的蝶形运算 
				for(k=0;k<distance;k++)//第k次蝶形运算  
					float theta=(float)(-2*PI*k/N);//旋转因子,WN(k)
					tempr1=Array[N*i+k];	
					tempi1=IArray[N*i+k];	
					tempr2=Math.cos(theta)*Array[N*i+k+distance]-Math.sin(theta)*IArray[N*i+k+distance];    //CB的实数				
					tempi2=Math.sin(theta)*Array[N*i+k+distance]+Math.cos(theta)*IArray[N*i+k+distance];    //CB的复数
					Array[N*i+k]=(float)(tempr1+tempr2);//(A+CB)的实数	
					IArray[N*i+k]=(float)(tempi1+tempi2);//(A+CB)的复数
					Array[N*i+k+distance]=(float)(tempr1-tempr2);//(A-CB)的实数
					IArray[N*i+k+distance]=(float)(tempi1-tempi2);//(A-CB)的复数
				 
			
			N=N*2;		//下一级蝶形运算中循环分母N*2
			distance*=2; //下一级蝶形运算两节点间的距离增加
			group/=2;  //下一级蝶形运算的组数减半
		
//
		//输出运算结果
		int t1=0;  //数组输出的循环次数
		System.out.println("");
		System.out.println("FFT运算结果如下：");
		while(t1<arraylength)
			 System.out.printf("X[%d]= %f\\t+j(%f)\\n",t1,Array[t1],IArray[t1]);
			 t1++;
		
	  


//阿汪先生的博客.ws

程序效果

1、复数输入

2、输入实部时，输入任一非数字字符（串），自动结束录入过程，剩余序列位自动补零；

在Matlab中验算

x=[1 2 3 0];
X=fft(x,4);
X

一些说明：

1、可复数输入，输入位数不足2^N次，自动补位；

2、输入实部时，输入任一非数字字符，自动结束录入过程，剩余序列位自动补零；

3、倒位运算采用了雷德算法；

4、计算机运算时位数的限制导致上述的两段程序结果的误差，问题不大；

5、修改因子W_N^K中的N可算任意N点的FFT运算，本程序中的N为(2 ^M)，具有一定局限性。

//只要在输入时单独赋值，即可实现序列任意N点（可奇数点）的FFT运算。