FFT---快速傅里叶变换

 

 最近在学习dsp数字信号处理算法,FFT是避免不了要学习的,本文将用实例理解FFT的实现。

FFT是DFT的快速计算方法,所以先简单的复习一下DFT的基本原理。

  DFT:

 

式中,,称为蝶形因子。上面的式子是N点的DFT,需要做的乘法是N^2,FFT算法是利用了蝶形因子WN的对称性和周期性,将乘法的运算量减少到N^2 /2。

 

蝶形因子的三个性质

        周期性:

        对称性: 和 

        可约性: 和 

 

 

 

FFT原理推导

 

 

 

 

 

 

 以N=8点的FFT计算为例。

注意:

1.第一层L=1时,N=2,;第二层L=2时,N=4;第三层L=3时,N=8.

2.红线表示“+”,蓝线表示“-”.

3.     i表示从左边出发原始的序列标号,j表示第j层。比如我标记绿色三角形,,表示x[1]的第二层。

4.最左边的x是小写的,右边的X是大写的.

 

 

 

 

 

 实例-----以计算X[1]为例

DFT的算法:每次算N个点

FFT的算法:每次算N/2个点

第一层:

 

 第二层:

 

 第三层:

 

然后通过递归迭代,把第二层带入第三层,把第一层带入第二层,最终算的X[1]。相当于把问题分解成了众多子问题,用算法描述就是T[N]=2*T[N/2]+cN,它的时间复杂度下降了,从而实现FFT的快速计算。

 

C语言实现:

复制代码
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#define PI 3.1415926535

void kfft(pr, pi, n, k, fr, fi) 
int n, k;
double pr[], pi[], fr[], fi[];
{
    int it, m, is, i, j, nv, l0;
    double p, q, s, vr, vi, poddr, poddi;
    for (it = 0; it <= n - 1; it++) //将pr[0]和pi[0]循环赋值给fr[]和fi[]
    {
        m = it;
        is = 0;
        for (i = 0; i <= k - 1; i++)
        {
            j = m / 2;
            is = 2 * is + (m - 2 * j);
            m = j;
        }
        fr[it] = pr[is];
        fi[it] = pi[is];
    }
    pr[0] = 1.0;
    pi[0] = 0.0;
    p = 6.283185306 / (1.0 * n);
    pr[1] = cos(p); //将w=e^-j2pi/n用欧拉公式表示
    pi[1] = -sin(p);

    for (i = 2; i <= n - 1; i++) //计算pr[]
    {
        p = pr[i - 1] * pr[1];
        q = pi[i - 1] * pi[1];
        s = (pr[i - 1] + pi[i - 1]) * (pr[1] + pi[1]);
        pr[i] = p - q;
        pi[i] = s - p - q;
    }
    for (it = 0; it <= n - 2; it = it + 2)
    {
        vr = fr[it];
        vi = fi[it];
        fr[it] = vr + fr[it + 1];
        fi[it] = vi + fi[it + 1];
        fr[it + 1] = vr - fr[it + 1];
        fi[it + 1] = vi - fi[it + 1];
    }
    m = n / 2;
    nv = 2;
    for (l0 = k - 2; l0 >= 0; l0--) //蝴蝶操作
    {
        m = m / 2;
        nv = 2 * nv;
        for (it = 0; it <= (m - 1) * nv; it = it + nv)
            for (j = 0; j <= (nv / 2) - 1; j++)
            {
                p = pr[m * j] * fr[it + j + nv / 2];
                q = pi[m * j] * fi[it + j + nv / 2];
                s = pr[m * j] + pi[m * j];
                s = s * (fr[it + j + nv / 2] + fi[it + j + nv / 2]);
                poddr = p - q;
                poddi = s - p - q;
                fr[it + j + nv / 2] = fr[it + j] - poddr;
                fi[it + j + nv / 2] = fi[it + j] - poddi;
                fr[it + j] = fr[it + j] + poddr;
                fi[it + j] = fi[it + j] + poddi;
            }
    }
    for (i = 0; i <= n - 1; i++)
    {
        pr[i] = sqrt(fr[i] * fr[i] + fi[i] * fi[i]); //幅值计算
    }
    return;
}

int main()
{
    int i, j;
    int N = 16;
    double pr[16] = {
        0.537667139546100, 1.83388501459509, -2.25884686100365, 0.862173320368121, 0.318765239858981, -1.30768829630527, -0.433592022305684, 0.342624466538650, 3.57839693972576, 2.76943702988488, -1.34988694015652, 3.03492346633185, 0.725404224946106, -0.0630548731896562, 0.714742903826096, -0.204966058299775};
    double pi[16] = {0.0}, fr[N], fi[N], t[N];

    // pr表示采样数据的实部,pi表示采样数据的虚部,fr表示离散傅里叶变换结果的实部,fi表示离散傅里叶变换结果的虚部。
    /*
    for (i=0; i<N; i++)  //生成输入信号
    {
        pr[i]=0.6*sin(2*PI*500*i)+0.6*sin(2*PI*50*i); pi[i]=0.0; //依次将值赋值给pr,pi赋值为0
    }
    */

    kfft(pr, pi, N,4, fr, fi, 0, 1); //调用FFT函数
    for (i = 0; i < N; i++)
    {
//        printf("%d\t%lf\n", i, pr[i]); //输出结果
        printf("%lf+(%lf)j\n",fr[i],fi[i]);
    }
    
}
复制代码

用matlab验证C代码

1
2
x = [0.537667139546100, 1.83388501459509, -2.25884686100365, 0.862173320368121, 0.318765239858981, -1.30768829630527, -0.433592022305684, 0.342624466538650, 3.57839693972576, 2.76943702988488, -1.34988694015652, 3.03492346633185, 0.725404224946106, -0.0630548731896562, 0.714742903826096, -0.204966058299775];
y = fft(x)

发现结果一致。

 

posted @   傅红雪a  阅读(1250)  评论(0编辑  收藏  举报
相关博文:
阅读排行:
· DeepSeek “源神”启动!「GitHub 热点速览」
· 微软正式发布.NET 10 Preview 1:开启下一代开发框架新篇章
· C# 集成 DeepSeek 模型实现 AI 私有化(本地部署与 API 调用教程)
· DeepSeek R1 简明指南:架构、训练、本地部署及硬件要求
· NetPad:一个.NET开源、跨平台的C#编辑器
Live2D
欢迎阅读『FFT---快速傅里叶变换』
点击右上角即可分享
微信分享提示