有用过Goertzel算法进行DFT做DTMF解码的吗？

xyzhang_yjs08

RT！！
    最近实验室里面要做DTMF解码，查阅了一些资料发现网上很多文章都说要用Gertzel算法进行DFT，然后根据其幅度谱来判断是否有某一频率的谱线。而且我认为这个算法本身应该是没有问题的，因为印象中好像PIC的库函数里面有这么一个解码的程序（不过我没有源代码）。
    但是实际上自己用这个变换做了分析之后发现，结果根本不是那么回事——计算出来的谱线根本就是乱的，根本无法区别。
    发个处理程序吧，希望有达人看到后可以指出问题的所在，程序是软件模拟的，数据是实际的采样值。
采样频率：8192Hz
DFT点数：206
采样点数：206
计算的频率应该是最前面的两个频率的合成，即：312.7和346.6那两个频率
程序比较简单就没加注释，呵呵


#include "stdio.h"
#include "math.h"
#define PI 3.1415926

short SourceData1[] =
{       
801,893,970,1012,1017,978,906,811,621,437,
289,65,-90,-180,-409,-512,-524,-513,-428,-187,
-113,-2,200,370,517,701,825,895,923,929,
885,812,660,493,347,173,-10,-114,-185,-398,
-447,-462,-436,-201,-170,-92,2,181,334,464,
589,687,743,804,800,698,623,503,389,293,
116,-18,-101,-155,-187,-196,-198,-185,-164,-114,
-42,68,169,301,369,441,493,515,515,507,
475,425,360,296,182,102,0,-36,-83,-102,
-105,-104,-91,-56,-26,7,86,142,195,284,
297,301,308,308,298,290,219,189,164,130,
108,98,77,80,73,80,90,92,101,110,
114,109,104,96,82,67,55,6,0,-4,
-4,2,51,75,109,152,184,223,291,305,
314,314,309,297,226,173,115,57,-25,-78,
-117,-145,-160,-166,-161,-148,-109,-63,-6,107,
206,311,378,437,480,493,484,460,396,321,
217,106,-28,-107,-179,-381,-434,-508,-509,-442,
-391,-182,-114,-25,132,296,400,495,603,664,
674,652,588,485,375,229,71,-76,-174,-412,
-540,-600,-629,-628,-580,-515,-374,
};

double frequence[16] =
{
    312.6,346.7,384.6,426.6,473.2,524.8,582.1,645.7,
    716.1,794.3,881.0,977.2,1083.9,1202.3,1333.5,1479.1
};

int main(void)
{
        int i,j;
        double dftres[32];
        double dftpow[16];
        double temp;
//DFT结果的第0点
        for(i=0;i<32;i+=2)
            dftres = (double)SourceData1[0];
//DFT中间结果的第1点
        for(i=0;i<32;i+=2)
            dftres[i+1] = 2.0*cos(2*PI*frequence[i/2]/8192.0)*dftres+(double)SourceData1[1];
        for(i=2;i<206;i++)
        {
                for(j=0;j<32;j+=2)
                {
                        temp = dftres[j];
                        dftres[j] = dftres[j+1];
                        dftres[j+1] = 2.0*cos(2*PI*frequence[j/2]/8192.0)*dftres[j]-temp+(double)SourceData1;
                }
        }
        for(j=0;j<16;j++)
                dftpow[j] = dftres[j]*dftres[j]+dftres[2*j+1]*dftres[2*j+1]
                                        -2.0*cos(2*PI*frequence[j]/8192.0)*dftres[j]*dftres[2*j+1];
        for(j=0;j<16;j++)
                printf("%16.2f\n",dftpow[j]);
        printf("\n");
        return 0;
}

comway

嘛

ldhavr

等待结果