| 网站首页 | 业界新闻 | 小组 | 威客 | 人才 | 下载频道 | 博客 | 代码贴 | 在线编程 | 编程论坛
欢迎加入我们,一同切磋技术
用户名:   
 
密 码:  
共有 1701 人关注过本帖
标题:用matlab自写fft函数,为什么和mat的fft结果有些差别?(附图,附代码 )
只看楼主 加入收藏
cacar
Rank: 1
等 级:新手上路
帖 子:1
专家分:0
注 册:2008-6-5
收藏
 问题点数:0 回复次数:0 
用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;
}
搜索更多相关主题的帖子: fft matlab 函数 附图 代码 
2008-06-05 10:11
快速回复:用matlab自写fft函数,为什么和mat的fft结果有些差别?(附图,附代码 ...
数据加载中...
 
   



关于我们 | 广告合作 | 编程中国 | 清除Cookies | TOP | 手机版

编程中国 版权所有,并保留所有权利。
Powered by Discuz, Processed in 0.019657 second(s), 9 queries.
Copyright©2004-2024, BCCN.NET, All Rights Reserved