FFT算法的验证正确性 [英] Verifying correctness of FFT algorithm
本文介绍了FFT算法的验证正确性的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！
问题描述
今天，我写了一个算法来计算从给定的分数排列的快速傅立叶变换重新presenting离散函数。现在，我想测试一下，看看它是否工作。我试着约十几个不同的输入集，他们似乎与我在网上找到的例子匹配。但是，对于我的最终测试，我给它的cos（I / 2）输入，与我从0到31，我已经得到了3种不同的结果，在此基础上求解器我用。我的解决方案似乎是最不准确的：
这是否表明我的算法有问题，或者是比较小的数据集的只是一个结果呢？
我的code是下面，在情况下，它可以帮助：
/ **
*片原始数组，先从开始，抓住每一个箭步元素。
*例如，切片（A，3，4，5）将返回从阵列A元件3，8，13，18和
* @参数阵列的阵列到被切
*参数启动的起始索引
* @参数满足newLength最后阵列的长度
* @参数元素跨步之间的间隔被选择
返回：输入数组的一个切片副本
* /
公共双[]切片（双[]数组，诠释开始，诠释满足newLength，INT步幅）{
双[] newArray =新的双[满足newLength]
诠释计数= 0;
的for（int i =启动; COUNT＆LT;满足newLength和放大器;＆安培; I＆LT; array.length; I + =步幅）{
newArray [统计++] =阵列[I]
}
返回newArray;
}
/ **
*计算快速傅立叶变换给定的功能。的参数与计算出的值来更新
*要忽略所有的虚输出，留下想象空
* @参数真实的重新$ P $阵列psenting的离散时间函数的实部
* @参数假想重新$ P $阵列psenting的离散时间函数的虚部
* pre：如果虚不空，这两个数组必须是相同的长度，它必须是2的幂
* /
公共无效FFT（双[]实数，双[]虚构的）抛出：IllegalArgumentException - {
如果（真正的== NULL）{
抛出新的NullPointerException异常（真正的数组不能为空）;
}
INT N = real.length;
//确保的长度是2的幂
如果（（将Math.log（N）/将Math.log（2））％1！= 0）{
抛出新抛出：IllegalArgumentException（数组的长度必须是2的幂）;
}
如果（虚= NULL和放大器;！＆安培;！imaginary.length = N）{
抛出新抛出：IllegalArgumentException（两个数组必须是相同的长度）;
}
如果（N == 1）{
返回;
}
双[] even_re =切片（真实，0，N / 2，2）;
双[] odd_re =片（真正的，1，N / 2，2）;
双[] even_im = NULL;
双[] odd_im = NULL;
如果（虚！= NULL）{
even_im =切片（虚，0，N / 2，2）;
odd_im =切片（虚，1，N / 2，2）;
}
FFT（even_re，even_im）;
FFT（odd_re，odd_im）;
// F [k]的实= [K] +假想[k]的
//偶奇
// F [K] = E [K] + O [K] * E ^（ - I * 2 * PI * K / N）
// F [K + N / 2] = E [K] - Ø[K] * E ^（ - I * 2 * PI * K / N）
//拆分复杂的阵列到阵列组件：
// E [K] = ER [K] + I * EI [K]
// -O [K] =或[K] + I * OI [K]
// E ^ IX = COS（X）+ I *的sin（x）
//令x = -2 * PI * K / N
// F [k]的呃= [K] + I *的ei [K] +（或[K] + I * OI [k]的）（COS（X）+ I *的sin（x））
// = ER [K] + I * EI [K] +或[K] COS（X）+ I *或[K]的sin（x）+ I * OI [K] COS（X） - OI [K]的sin（x）
// =（ER [K] +或[K]为cos（x）的 - OI [k]的的sin（x））+ I *（EI [K] +或[K]的sin（x）+ OI [k]的余弦（x）的）
// {实} {}虚
// F〔K + N / 2] =（ER [K] - 或[k]的余弦（x）的+ OI [k]的的sin（x））+ I *（EI [K] - 或[k]的罪（ x）的 - OI [k]的余弦（x）的）
// {实} {}虚
//忽略所有的虚部（OI = 0）：
// F [k]的呃= [K] +或[k]的余弦（x）的
// F〔K + N / 2] =呃[K] - 或[k]的余弦（x）的
为（中间体K = 0; K＆所述N / 2 ++ k）的{
双T = odd_re [K] * Math.cos（-2 * Math.PI * K / N）;
真实[K] = even_re [K] +吨;
真正的[K + N / 2] = even_re [K] - 吨;
如果（虚！= NULL）{
T = odd_im [K] * Math.sin（-2 * Math.PI * K / N）;
真正的[K] - = T;
真正的[K + N / 2] + = T;
双T1 = odd_re [K] * Math.sin（-2 * Math.PI * K / N）;
双T2 = odd_im [K] * Math.cos（-2 * Math.PI * K / N）;
假想[K] = even_im [K] + T1 + T2;
虚[K + N / 2] = even_im [K] - T1 - T2;
}
}
}
解决方案
1.Validation
- 看这里：慢DFT和IDFT
- 在到底是我的缓慢实施DFT和的iDFT的
- 测试和正确的我也过去 使用它快速实现验证
2.您code
- 停止递归是错误的（你忘了设置返回元素）
-
我的是这样的：
如果（N＆LT; = 1）{如果（N == 1）{DST [0] = SRC [0] * 2.0; DST [1] = SRC [1] * 2.0; } 返回; } -
所以当你的ñ== 1设置输出元素重新= 2.0 *真正的[0]，林= 2.0 *虚[0]返回之前。
- 还我有点在复杂的数学（T，T1，T2）和懒惰输给了分析
- ，但只是要确定这里是我的快速实施
- 需要从类层次太多的事情
- ，所以它不会是另一个需要你了，然后视觉比较你的code
我的快速实现（CC意味着复杂的输出和输入）：
// ------------------------------------ ---------------------------------------
无效变换:: DFFTcc（双* DST，双* SRC，INT N）
{
如果（正将N）的init（n）的;
如果（正＆其中; = 1）{如果（N == 1）{DST [0] = SRC [0] * 2.0; DST [1] = SRC [1] * 2.0; } 返回; }
INT I，J，N = N＆GT;＆GT; 1，Q，DQ = + N / N，MQ = N-1;
//重新排序偶数，奇数（buterfly）
为（J = 0，I = 0; I＆所述; N + N）{DST [J] = SRC [i]于;我++; J ++; DST [J] SRC = [I] I + = 3; J ++; }
为（ⅰ= 2; I＆所述; N + N）{DST [J] = SRC [i]于;我++; J ++; DST [J] SRC = [I] I + = 3; J ++; }
//递归
DFFTcc（SRC，DST，N2）; // 甚至
DFFTcc（SRC + N，DST + N，N2）; // 奇
//重新排序和体重恢复（buterfly）
双A0，A1，B0，B1，A，B;
为（Q = 0，I = 0，J =正;我n种; I + = 2，J + = 2，Q =（Q + DQ）及MQ）
{
A0 = SRC [J]。 A1 = + _ COS [Q]
B0 = SRC [J + 1]; B1 = + _罪[Q]
一个=（A0 * A1） - （B0·B1）;
B =（A0 * B1）+（A1 * B0）;
A0 = SRC [I] A1 =一个;
B0 = SRC [I + 1]; B1 = B;
DST [I] =（A0 + A1）* 0.5;
DST [I + 1] =（B0 + B1）* 0.5;
DST [J] =（A0-A1）* 0.5;
DST [J + 1] =（B0-B1）* 0.5;
}
}
// ------------------------------------------------ ---------------------------
- 在DST []和src []不重叠的！
- ，所以你不能改变阵列本身
- _cos和_sin是precomputed COS和罪恶值的表（由init（）计算功能
-
是这样的：
双A，DA; INT I; DA = 2.0 * M_PI /双（N）; 用于：（a = 0.0，I = 0; I＆所述N;我+ +，一个+ =哒）{_cos [I] = COS（一）; _sin [I] =罪（一）; } -
N为2（数据集的补零大小）功率（从INIT（n）的最后调用N）
只是要完成这里是我的复杂，复杂的慢版：
// ------------------------------------ ---------------------------------------
无效变换:: DFTcc（双* DST，双* SRC，INT N）
{
INT I，J;
双重的a，b，A0，A1，_n，B0，B1，Q，QQ，DQ;
DQ = + 2.0 * M_PI /双（N）; _n = 2.0 /双（N）;
为（Q = 0.0，J = 0; J＆n种; J ++，Q + = DQ）
{
A = 0.0; B = 0.0;
对于（QQ = 0.0，I = 0;我n种;我++，QQ + = Q）
{
A0 = SRC [I + I]。 A1 = + COS（QQ）;
B0 = SRC [I + I + 1]; B1 = +罪（QQ）;
一个+ =（A0 * A1） - （B0·B1）;
B + =（A0 * B1）+（A1 * B0）;
}
DST [J + J] = A * _n;
DST [J + J + 1] = B * _n;
}
}
// ------------------------------------------------ ---------------------------
Today I wrote an algorithm to compute the Fast Fourier Transform from a given array of points representing a discrete function. Now I'm trying to test it to see if it is working. I've tried about a dozen different input sets, and they seem to match up with examples I've found online. However, for my final test, I gave it the input of cos(i / 2), with i from 0 to 31, and I've gotten 3 different results based on which solver I use. My solution seems to be the least accurate:
Does this indicate a problem with my algorithm, or is it simply a result of the relatively small data set?
My code is below, in case it helps:
/**
* Slices the original array, starting with start, grabbing every stride elements.
* For example, slice(A, 3, 4, 5) would return elements 3, 8, 13, and 18 from array A.
* @param array The array to be sliced
* @param start The starting index
* @param newLength The length of the final array
* @param stride The spacing between elements to be selected
* @return A sliced copy of the input array
*/
public double[] slice(double[] array, int start, int newLength, int stride) {
double[] newArray = new double[newLength];
int count = 0;
for (int i = start; count < newLength && i < array.length; i += stride) {
newArray[count++] = array[i];
}
return newArray;
}
/**
* Calculates the fast fourier transform of the given function. The parameters are updated with the calculated values
* To ignore all imaginary output, leave imaginary null
* @param real An array representing the real part of a discrete-time function
* @param imaginary An array representing the imaginary part of a discrete-time function
* Pre: If imaginary is not null, the two arrays must be the same length, which must be a power of 2
*/
public void fft(double[] real, double[] imaginary) throws IllegalArgumentException {
if (real == null) {
throw new NullPointerException("Real array cannot be null");
}
int N = real.length;
// Make sure the length is a power of 2
if ((Math.log(N) / Math.log(2)) % 1 != 0) {
throw new IllegalArgumentException("The array length must be a power of 2");
}
if (imaginary != null && imaginary.length != N) {
throw new IllegalArgumentException("The two arrays must be the same length");
}
if (N == 1) {
return;
}
double[] even_re = slice(real, 0, N/2, 2);
double[] odd_re = slice(real, 1, N/2, 2);
double[] even_im = null;
double[] odd_im = null;
if (imaginary != null) {
even_im = slice(imaginary, 0, N/2, 2);
odd_im = slice(imaginary, 1, N/2, 2);
}
fft(even_re, even_im);
fft(odd_re, odd_im);
// F[k] = real[k] + imaginary[k]
// even odd
// F[k] = E[k] + O[k] * e^(-i*2*pi*k/N)
// F[k + N/2] = E[k] - O[k] * e^(-i*2*pi*k/N)
// Split complex arrays into component arrays:
// E[k] = er[k] + i*ei[k]
// O[k] = or[k] + i*oi[k]
// e^ix = cos(x) + i*sin(x)
// Let x = -2*pi*k/N
// F[k] = er[k] + i*ei[k] + (or[k] + i*oi[k])(cos(x) + i*sin(x))
// = er[k] + i*ei[k] + or[k]cos(x) + i*or[k]sin(x) + i*oi[k]cos(x) - oi[k]sin(x)
// = (er[k] + or[k]cos(x) - oi[k]sin(x)) + i*(ei[k] + or[k]sin(x) + oi[k]cos(x))
// { real } { imaginary }
// F[k + N/2] = (er[k] - or[k]cos(x) + oi[k]sin(x)) + i*(ei[k] - or[k]sin(x) - oi[k]cos(x))
// { real } { imaginary }
// Ignoring all imaginary parts (oi = 0):
// F[k] = er[k] + or[k]cos(x)
// F[k + N/2] = er[k] - or[k]cos(x)
for (int k = 0; k < N/2; ++k) {
double t = odd_re[k] * Math.cos(-2 * Math.PI * k/N);
real[k] = even_re[k] + t;
real[k + N/2] = even_re[k] - t;
if (imaginary != null) {
t = odd_im[k] * Math.sin(-2 * Math.PI * k/N);
real[k] -= t;
real[k + N/2] += t;
double t1 = odd_re[k] * Math.sin(-2 * Math.PI * k/N);
double t2 = odd_im[k] * Math.cos(-2 * Math.PI * k/N);
imaginary[k] = even_im[k] + t1 + t2;
imaginary[k + N/2] = even_im[k] - t1 - t2;
}
}
}
解决方案
1.Validation
- look here: slow DFT,iDFT
- at the end is mine slow implementation of DFT and iDFT
- tested and correct I also use it for fast implementations validation in the past
2.Your code
- stop recursion is wrong (you forget to set the return element)
mine looks like this:
if (n<=1) { if (n==1) { dst[0]=src[0]*2.0; dst[1]=src[1]*2.0; } return; }so when your N==1 set the output element to Re=2.0*real[0], Im=2.0*imaginary[0] before return.
- also I am a bit lost in your complex math (t,t1,t2) and to lazy to analyze
- but just to be sure here is mine fast implementation
- need too much things from class hierarchy
- so it will not be of another use for you then visual comparison to your code
My Fast implementation (cc means complex output and input):
//---------------------------------------------------------------------------
void transform::DFFTcc(double *dst,double *src,int n)
{
if (n>N) init(n);
if (n<=1) { if (n==1) { dst[0]=src[0]*2.0; dst[1]=src[1]*2.0; } return; }
int i,j,n2=n>>1,q,dq=+N/n,mq=N-1;
// reorder even,odd (buterfly)
for (j=0,i=0;i<n+n;) { dst[j]=src[i]; i++; j++; dst[j]=src[i]; i+=3; j++; }
for ( i=2;i<n+n;) { dst[j]=src[i]; i++; j++; dst[j]=src[i]; i+=3; j++; }
// recursion
DFFTcc(src ,dst ,n2); // even
DFFTcc(src+n,dst+n,n2); // odd
// reorder and weight back (buterfly)
double a0,a1,b0,b1,a,b;
for (q=0,i=0,j=n;i<n;i+=2,j+=2,q=(q+dq)&mq)
{
a0=src[j ]; a1=+_cos[q];
b0=src[j+1]; b1=+_sin[q];
a=(a0*a1)-(b0*b1);
b=(a0*b1)+(a1*b0);
a0=src[i ]; a1=a;
b0=src[i+1]; b1=b;
dst[i ]=(a0+a1)*0.5;
dst[i+1]=(b0+b1)*0.5;
dst[j ]=(a0-a1)*0.5;
dst[j+1]=(b0-b1)*0.5;
}
}
//---------------------------------------------------------------------------
- dst[] and src[] are not overlapping !!!
- so you cannot transform array to itself
- _cos and _sin are precomputed tables of cos and sin values (computed by init() function
like this:
double a,da; int i; da=2.0*M_PI/double(N); for (a=0.0,i=0;i<N;i++,a+=da) { _cos[i]=cos(a); _sin[i]=sin(a); }N is power of 2 (zero padded size of data set) (last n from init(n) call)
Just to be complete here is mine complex to complex slow version:
//---------------------------------------------------------------------------
void transform::DFTcc(double *dst,double *src,int n)
{
int i,j;
double a,b,a0,a1,_n,b0,b1,q,qq,dq;
dq=+2.0*M_PI/double(n); _n=2.0/double(n);
for (q=0.0,j=0;j<n;j++,q+=dq)
{
a=0.0; b=0.0;
for (qq=0.0,i=0;i<n;i++,qq+=q)
{
a0=src[i+i ]; a1=+cos(qq);
b0=src[i+i+1]; b1=+sin(qq);
a+=(a0*a1)-(b0*b1);
b+=(a0*b1)+(a1*b0);
}
dst[j+j ]=a*_n;
dst[j+j+1]=b*_n;
}
}
//---------------------------------------------------------------------------
这篇关于FFT算法的验证正确性的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
查看全文
