打印

fft采样的一点疑惑

[复制链接]
2230|2
手机看帖
扫描二维码
随时随地手机跟帖
跳转到指定楼层
楼主
sddds|  楼主 | 2011-3-3 11:42 | 只看该作者 回帖奖励 |倒序浏览 |阅读模式
沙发
sddds|  楼主 | 2011-3-3 17:47 | 只看该作者
查了一下资料 应该是任意起始点都可以
还有一点疑惑 FFT计算时正弦表余弦表怎么来的

使用特权

评论回复
板凳
aresc| | 2011-3-3 23:36 | 只看该作者
void RADIX_2_DIT_FFT(double *fftBuffer, long fftFrameSize, long sign)
/*
        FFT routine, (C)1996 S.M.Bernsee. Sign = -1 is FFT, 1 is iFFT (inverse)
        Fills fftBuffer[0...2*fftFrameSize-1] with the Fourier transform of the
        time domain data in fftBuffer[0...2*fftFrameSize-1]. The FFT array takes
        and returns the cosine and sine parts in an interleaved manner, ie.
        fftBuffer[0] = cosPart[0], fftBuffer[1] = sinPart[0], asf. fftFrameSize
        must be a power of 2. It expects a complex input signal (see footnote 2),
        ie. when working with 'common' audio signals our input signal has to be
        passed as {in[0],0.,in[1],0.,in[2],0.,...} asf. In that case, the transform
        of the frequencies of interest is in fftBuffer[0...fftFrameSize].
        2-radix DIT(Decimation In Time)-FFT algorithm.
*/
{
        double wr, wi, arg, *p1, *p2, temp;
        double tr, ti, ur, ui, *p1r, *p1i, *p2r, *p2i;
        long i, bitm, j, le, le2, k;

        // re-order x(n) or X[n] in bit reverse order.
        // both x(n) and X[n] are complex value, a complex value occupies contineous
        // two memory, the first is real part, and the second is imaginary part.
        for (i = 2; i < 2*fftFrameSize-2; i += 2) {
                for (bitm = 2, j = 0; bitm < 2*fftFrameSize; bitm <<= 1) {
                        if (i & bitm) j++;
                        j <<= 1;
                }

                // swap two point according bit reverse order
                if (i < j) {
                        p1 = fftBuffer+i; p2 = fftBuffer+j;
                        temp = *p1; *(p1++) = *p2;
                        *(p2++) = temp; temp = *p1;
                        *p1 = *p2; *p2 = temp;
                }
        }

        // After re-order in bit reverse rule, then the output will be in normal order.
        for (k = 0, le = 2; k < (long)(log(fftFrameSize)/log(2.)+.5); k++) {
                le <<= 1;
                le2 = le>>1;
                ur = 1.0;                                // cos0
                ui = 0.0;                                // sin0
                arg = M_PI / (le2>>1);        // 2*PI/(2^r)
                wr = cos(arg);                        // wr = cos(arg)
                wi = sign*sin(arg);                // wi = sin(sign*arg), if forward FFT, then sign = -1, otherwise sign = 1
                for (j = 0; j < le2; j += 2) {
                        p1r = fftBuffer+j; p1i = p1r+1;                // the first point of one butterfly calculation: P1
                        p2r = p1r+le2; p2i = p2r+1;                        // the second point of one butterfly calculation: P2
                        for (i = j; i < 2*fftFrameSize; i += le) {
                                tr = *p2r * ur - *p2i * ui;                // real part of: P2 * WN^r
                                ti = *p2r * ui + *p2i * ur;                // imag part of: P2 * WN^r
                                *p2r = *p1r - tr; *p2i = *p1i - ti;
                                *p1r += tr; *p1i += ti;
                                p1r += le; p1i += le;
                                p2r += le; p2i += le;
                        }
                        tr = ur*wr - ui*wi;                // cos(A+arg) = cosAcos(arg) - sinAsin(arg)
                        ui = ur*wi + ui*wr;                // sin(A+arg) = sinAcos(arg) + cosAsin(arg)
                        ur = tr;
                }
        }
}

使用特权

评论回复
发新帖 我要提问
您需要登录后才可以回帖 登录 | 注册

本版积分规则

4

主题

21

帖子

1

粉丝