Java implementation Fast Fourier Transform FFT (Base-2 FFT for Time Extraction, Inverse Calculation, Butterfly Operation)
import java.util.Scanner; import java.math.*; //log function class Log{ static public double log(double value, double base) { return Math.log(value) / Math.log(base); } } //power function 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){ //init varaibles final double PI=3.141; int N1=0; //series length int M=0; //series level 2^M int arraylength; // int i=0; //input array loop count int t=0; //output array loop count /////////////////////////////////////////////////////////////// //check Scanner sc=new Scanner (System.in); while(N1<2){ System.out.println("check your series length!"); 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++; } } }
Program effect
1.Plural input
2. When entering the real part, enter any non-numeric character (string) to automatically end the entry process and the remaining sequence bits are automatically zero-added;
Checking in Matlab
x=[1 2 3 0]; X=fft(x,4); X
Some notes:
1. Plural input is possible. If the number of input digits is less than 2 ^ N times, it will be filled automatically.
2. When entering the real part, enter any non-numeric character to automatically end the entry process and zero-padded the remaining sequence bits automatically;
3. The inversion operation uses the Reid algorithm;
4. The limitation of the number of digits in the computer operation results in the error of the above two programs, and the problem is not large;
5. Modify the N in the factor W_N ^ K to calculate the FFT operation of any N points. The N in this program is (2 ^ M), which has certain limitations.
// As long as the input is individually assigned, the FFT operation of any N points (may be odd points) in the sequence can be realized.
Orignal link:https://blog.csdn.net/qq_43499622/article/details/102759545