用matlab自写fft函数,为什么和mat的fft结果有些差别?(附图,附代码 )
是个作业,照cpp的fft代码照搬了一个,做出来和mat自己的fft不太一样,看不出来是哪的原因导致的差别,望指教~ 说明:对一个长度为128的序列(1:10=1,其他为0)进行fft,左边是mat的结果,右边是我的。。。
" border="0" />
我写的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; }