FFT.rar_fft vc_vc fft

// // 在本文件中, // (1) 定义了以直角坐标表示的复数和以极坐标表示的复数; // (2) 给出快速付里叶变换所需的子程序. // #include <stdio.h> #include <math.h> // 结构体 --- 以直角坐标表示的复数. struct complex { double re; // 复数的衧实部. double im; // 复数的虚部. }; // 结构体 --- 以极坐标表示的复数. struct rad { double amp; // 复数的幅值. double phase; // 复数的相角(以弧度表示). }; const double Pi = 3.14159265; complex wArray[2048]; // 计算 complex 型复数 c1 (被加数) 与 c2 (加数) 之和. 和数在函数的返回值中. complex cadd(complex c1, complex c2) { complex temp; temp.re = c1.re + c2.re; temp.im = c1.im + c2.im; return temp; } // 计算 complex 型复数 c1 (被减数) 与 c2 (减数) 之差. 差数在函数的返回值中. complex csub(complex c1, complex c2) { complex temp; temp.re = c1.re - c2.re; temp.im = c1.im - c2.im; return temp; } // 计算 complex 型复数 c1 (被乘数) 与 c2 (乘积) 之乘积. 乘积在函数的返回值中. complex cmul(complex c1, complex c2) { complex temp; temp.re = c1.re * c2.re - c1.im * c2.im; temp.im = c1.re * c2.im + c1.im * c2.re; return temp; } // 计算 complex 型复数 c1 (被除数) 与 c2 (除数) 之商. 商在函数的返回值中. complex cdiv(complex c1, complex c2) { double k; complex temp1, temp2; temp2.re = c2.re; temp2.im = -c2.im; temp1 = cmul(c1, temp2); k = c2.re * c2.re + c2.im * c2.im; temp1.re = temp1.re / k; temp1.im = temp1.im / k; return temp1; } // 求 complex 型复数 c 的共轭值. 函数 cconju 返回该共轭值. complex cconju(complex c) { complex temp; temp.re = c.re; temp.im = -c.im; return temp; } // 计算 rad 型复数 r1 (被加数) 与 r2 (加数) 之和. 和数在函数的返回值中. rad radd(rad r1, rad r2) { double x, y; rad temp; x = r1.amp * cos(r1.phase) + r2.amp * cos(r2.phase); y = r1.amp * sin(r1.phase) + r2.amp * sin(r2.phase); temp.amp = sqrt(x * x + y * y); temp.phase = atan(y / x); return temp; } // 计算 rad 型复数 r1 (被减数) 与 r2 (减数) 之差. 差数在函数的返回值中. rad rsub(rad r1, rad r2) { double x, y; rad temp; x = r1.amp * cos(r1.phase) - r2.amp * cos(r2.phase); y = r1.amp * sin(r1.phase) - r2.amp * sin(r2.phase); temp.amp = sqrt(x * x + y * y); temp.phase = atan(y / x); return temp; } // 计算 rad 型复数 r1 (被乘数) 与 r2 (乘数) 之乘积. 乘积在函数的返回值中. rad rmul(rad r1, rad r2) { rad temp; temp.amp = r1.amp * r2.amp; temp.phase = r1.phase + r2.phase; return temp; } // 计算 rad 型复数 r1 (被除数) 与 r2 (除数) 之乘积. 商在函数的返回值中. rad rdiv(rad r1, rad r2) { rad temp; temp.amp = r1.amp / r2.amp; temp.phase = r1.phase - r2.phase; return temp; } // 求 rad 型复数 c 的共轭值. 函数 cconju 返回该共轭值. rad rcomju(rad r) { rad temp; temp.amp = r.amp; temp.phase = -r.phase; return temp; } // 将 complex 型复数 c 转换为 rad 型复数(在函数的返回值中). rad r_c(complex c) { rad temp; temp.amp = sqrt(c.re * c.re + c.im * c.im); temp.phase = atan(c.im / c.re); return temp; } // 将 rad 型复数 r 转换为 compkex 型复数(在函数的返回值中). complex c_r(rad r) { complex c; c.re = r.amp * cos(r.phase); c.im = r.amp * sin(r.phase); return c; } // 检查快速付里叶变换点数. 点数是 2 的整数幂时, 返回值是迭代次数, 否则, 返回值是 0.�. int checkNum(int num) { int iterNum = 0; do { if(num % 2==0) { num = num / 2; iterNum = iterNum + 1; if(num==2) { iterNum = iterNum + 1; return iterNum; } } else return 0; }while(1); } // 求迭代次数. num 是点数, 函数 iterNum 的返回值是迭代次数. int iterNum(int num) { int i; i = 0; do { if(num%2==0) { num = num / 2; i = i + 1; if(num==2) { i = i + 1; break; } } else return 0; }while(num!=2); return i; } // 求一个数的二进逆序值. iterNum --- 迭代次数; ord --- -整数. // 函数 bitReversed 的返回值是 ord 的二进逆序. int bitReversed(int iterNum, int ord) { int i, inc, rev; int intArray[15]; inc =1; for(i=0;i<=iterNum-1;i++) { intArray[iterNum-1-i] = inc; inc = inc * 2; } rev = 0; for(i=0;i<=iterNum-1;i++) { if(ord%2==1) rev = rev + intArray[i]; ord = ord / 2; } return rev; } // 计算旋转因子序列 W. number 是点数; FgCal 和 FgNeg 是标志; outArray[] 是计算所得的旋转因子 // 序列 W. void calW(int number, int FgCal, int FgNeg, complex outArray[]) { int i; double phase, temp, Pi; Pi = 3.14159265; if(FgCal==1) // FgCal = 1 时, 需要计算数组 outArray[] 的实部和虚部. { temp = 2 * Pi / number; // 等分圆周,点数为 number. outArray[0].re = 1.0; outArray[0].im = 0.0; phase = 0.0; for(i=1;i<=number-1;i++) { if(FgNeg==0) phase = phase - temp; // FgNeg = 0 时, 计算 exp(-phase). else phase = phase + temp; // FgNeg = 1 时, 计算 exp(+phase). outArray[i].re = cos(phase); outArray[i].im = sin(phase); } } if(FgCal==0) // FgCal = 0 时, 不需要计算数组 outArray[] 的实部和虚部. for(i=0;i<=number-1;i++) // 这时, 只将的虚部改变符号即可. outArray[i].im = -outArray[i].im; } // 双精度数组按二进逆序重排子程序. points 是点数; iterNum 是迭代次数; inOut[] 是双精度输入数组, // 运行这子程序后, inOut[] 成为按二进逆序重排的双精度数组. void reRange1( int points, int iterNum, double inOut[]) { int i, p; double temp; for(i=0;i<=points-1;i++) { p = bitReversed(iterNum, i); if(p>i) { temp = inOut[p]; inOut[p] = inOut[i]; inOut[i] = temp; } } } // 复数数组按二进逆序重排子程序. points 是点数; iterNum 是迭代次数; inOut[] 是复数输入数组, // 运行这子程序后, inOut[] 成为按二进逆序重排的复数数组. void reRange2( int points, int iterNum, complex inOut[]) { int i, p; complex temp; for(i=0;i<=points-1;i++) { p = bitReversed(iterNum, i); if(p>i) { temp.re = inOut[p].re; temp.im = inOut[p].im; inOut[p].re = inOut[i].re; inOut[p].im = inOut[i].im; inOut[i].re = temp.re; inOut[i].im = temp.im; } } } // 这个子程序将一个 complex 型复数 acom 变为 rad 型复数(在返回值中) ,其相位的绝对值不大于 Pi. rad complex_rad(complex acom) { rad arad; int fg0; fg0 = 0; arad.amp = sqrt((acom.re * acom.re) + (acom.im * acom.im)); if(fabs(acom.re) < 0.0000001) acom.re = 0; if(fabs(acom.im) < 0.0000001) acom.im = 0; if((acom.re == 0) && (acom.im > 0)) { arad.phase = Pi / 2; fg0 = 1; } else if((acom.re == 0) && (acom.im < 0)) { arad.phase = -Pi / 2; fg0 = 1; } else if((acom.re < 0) && (acom.im == 0)) { arad.phase = Pi; fg0 = 1; } if(fg0==0) { if(acom.re!=0) arad.phase = fabs(atan(acom.im / acom.re)); if((acom.re >= 0) && (acom.im >= 0)) arad.phase = arad.phase; if((acom.re < 0) && (acom.im > 0)) arad.phase = Pi - arad.phase; if((acom.re > 0) && (acom.im < 0)) arad.phase = -arad.phase; if((acom.re < 0) && (acom.im < 0)) arad.phase = Pi + arad.phase; } return arad; } /* FFT0A 是快速付里叶变换子程序,采用第一种时分算法： 输入序列按自然顺序排列,输出序列按倒码顺序排列. 在 FFT0A 子程序中， points 为点数； inOut（）为输入/输出同址序列； T 为时域序列采样间隔. 计算严格定义的 DFT 时,取 T = 1.0. 通过 FFT 计算近似的付里叶变换时, T 是所选的时域序列采样间隔 ; Flag01 为标志， Flag01 = 0 时，使用已经算好的旋转因子序列 W（）,不作任何改变 ; Flag01 = 1 时，计算 W（）； Flag02 为标志， Flag02 = 0 时，不重新计算 W（）序列，只将 W （）序列的虚部改变符号， Flag02 = 1 时，重新计算 W（）序列，这时，使用标志 Flag03； Flag03 为标志， Flag03 = 0 时，计算 W（-j phase），进行正变换， Flag03 = 1 时，计算 W（+j phase），进行反变换； Flag04 为标志， Flag04 = 0 时，只得到以直角座标表示的正变换序列（inOut( )）， Flag04 = 1 时，除了得到以直角座标表示的正变换序列（in