FFT算法的验证正确性 [英] Verifying correctness of FFT algorithm
问题描述
今天,我写了一个算法来计算从给定的分数排列的快速傅立叶变换重新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<满足newLength和放大器;&安培; I< 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< = 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>> 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屋!