PDA

查看完整版本 : [求助]用matlab自写fft函数,为什么和mat的fft结果有些差别?(附图,附代码 )


cacar
2008-06-05, 10:07
是个作业,照cpp的fft代码照搬了一个,做出来和mat自己的fft不太一样,看不出来是哪的原因导致的差别,望指教~
说明:对一个长度为128的序列(1:10=1,其他为0)进行fft,左边是mat的结果,右边是我的。。。
http://bloodfalcon.ik8.com/1.jpg

我写的mat代码:
function FD = exfft(TD,r)
%----------------------------------
% Name: exfft
% Param:
% TD -时域序列
% r -2的幂数,即迭代次数
% Return:
% FD -频域序列
% Desc: 对给定的序列进行快速傅立叶变换
% By:
%-----------------------------------
count = 2^r; % 计算变换点数

% 分别用来保存蝶形运算的输入和输出
X1 = zeros(1,2^r);
X2 = zeros(1,2^r);


% 计算加权系数
W = zeros(1,count/2);
for j=1:count/2
angle = j*3.1415926*2/count;
W(1,j)=cos(angle)-i*sin(angle);
end

% 采用蝶形算法进行快速傅立叶变换
X1=TD;
for k=1:r % 迭代次数
for j=1:(2^k) % 该次迭代有 2*k 个部分
bfSize = 2^(r-k); % 每部分的点数
for l=1:bfSize/2
p=(j-1)*bfSize;
X2(1,l+p) = X1(1,l+p) + X1(1,l+p+bfSize/2);
X2(1,l+p+bfSize/2) = (X1(1,l+p) - X1(1,l+p+bfSize/2))*W(1,l*(2^k));
end
end
% 将结果X2作为下一次迭代的输入
X=X1;
X1=X2;
X2=X;
end
% 重新排序
FD=bitrevorder(X2)
原先的cpp代码(这个代码应该没有问题):
/*************************************************************************
* 函数名称:
* FFT()
* 参数:
* complex<double> * TD - 指向时域数组的指针
* complex<double> * FD - 指向频域数组的指针
* r -2的幂数,即迭代次数
* 返回值:
* 无。
* 说明:
* 该函数用来实现快速付立叶变换。
************************************************************************/
VOID CDibImage::FFT(complex<double> * TD, complex<double> * FD, int r)
{
LONG count; // 付立叶变换点数
int i,j,k; // 循环变量
int bfsize,p;
double angle; // 角度
complex<double> *W,*X1,*X2,*X;

count = 1 << r; // 计算付立叶变换点数

// 分配运算所需存储器
W = new complex<double>[count / 2];
X1 = new complex<double>[count];
X2 = new complex<double>[count];

// 计算加权系数
for(i = 0; i < count / 2; i++)
{
angle = -i * PI * 2 / count;
W[i] = complex<double> (cos(angle), sin(angle));
}

// 将时域点写入X1
memcpy(X1, TD, sizeof(complex<double>) * count);

// 采用蝶形算法进行快速付立叶变换
for(k = 0; k < r; k++)
{
for(j = 0; j < 1 << k; j++)
{
bfsize = 1 << (r-k);
for(i = 0; i < bfsize / 2; i++)
{
p = j * bfsize;
X2[i + p] = X1[i + p] + X1[i + p + bfsize / 2];
X2[i + p + bfsize / 2] = (X1[i + p] - X1[i + p + bfsize / 2])
* W[i * (1<<k)];
}
}
X = X1;
X1 = X2;
X2 = X;
}

// 重新排序
for(j = 0; j < count; j++)
{
p = 0;
for(i = 0; i < r; i++)
{
if (j&(1<<i))
{
p+=1<<(r-i-1);
}
}
FD[j]=X1[p];
}

delete W;
delete X1;
delete X2;
}

cacar
2008-06-06, 21:32
下标错了,已解决

xu2006
2008-10-12, 13:33
楼主,你能说清楚一下你的问题吗?