FFT-C語言

時間  2019-11-10
標籤 fft c語言 简体版
原文   原文鏈接

http://blog.miskcoo.com/2015/04/polynomial-multiplication-and-fast-fourier-transform#i-21算法

1、完全理解傅里葉變換

　　快速傅里葉變換(Fast Fourier Transform)是離散傅里葉變換的一種快速算法，簡稱FFT，經過FFT能夠將一個信號從時域變換到頻域。編程

　　模擬信號通過A/D轉換變爲數字信號的過程稱爲採樣。爲保證採樣後信號的頻譜形狀不失真，採樣頻率必須大於信號中最高頻率成分的2倍，這稱之爲採樣定理。函數

假設採樣頻率爲fs，採樣點數爲N，那麼FFT結果就是一個N點的複數，每個點就對應着一個頻率點，某一點n(n從1開始)表示的頻率爲：fn=(n-1)*fs/N。spa

舉例說明：用1kHz的採樣頻率採樣128點，則FFT結果的128個數據即對應的頻率點分別是0，1k/128，2k/128，3k/128，…，127k/128 Hz。3d

這個頻率點的幅值爲：該點複數的模值除以N/2(n=1時是直流份量，其幅值是該點的模值除以N)。code

2、傅里葉變換的C語言編程　　

一、碼位倒序。

　　假設一個N點的輸入序列，那麼它的序號二進制數位數就是t=log2N.
　　碼位倒序要解決兩個問題：①將t位二進制數倒序;②將倒序後的兩個存儲單元進行交換。
　　若是輸入序列的天然順序號i用二進制數表示，例如若最大序號爲15，即用4位就可表示n3n2n1n0，則其倒序後j對應的二進制數就是n0n1n2n3，那麼怎樣才能實現倒序呢?利用C語言的移位功能!orm

/*變址計算，將x(n)碼位倒置*/

void Reverse(complex *IN_X, int n)
{
    int row = log(n) / log(2); //row爲FFT的N個輸入對應的最大列數
    complex temp;  //臨時交換變量
    unsigned short i = 0, j = 0, k = 0;
    double t;
    for (i = 0; i < row; i++)  //從第0個序號到第N-1個序號
    {
        k = i;   //當前第i個序號
        j = 0;   //存儲倒序後的序號，先初始化爲0
        t = row;  //共移位t次
        while ((t--) > 0)    //利用按位與以及循環實現碼位顛倒
        {
            j = j << 1;
            j |= (k & 1);    //j左移一位而後加上k的最低位
            k = k >> 1;      //k右移一位，次低位變爲最低位
                             /*j爲倒序，k爲正序，
                               從k右向左的值依次移位，
                               j向左依次填入對應的k的右向位*/
        }
        if (j > i)   //若是倒序後大於原序數，就將兩個存儲單元進行交換(判斷j>i是爲了防止重複交換
        {
            temp     = IN_X[i];
            IN_X[i] = IN_X[j];
            IN_X[j] = temp;
        }
    }
}


//複數加法的定義
complex add(complex a, complex b)  
{
    complex c;
    c.real = a.real + b.real;
    c.img = a.img + b.img;
    return c;
}

//複數乘法的定義
complex mul(complex a, complex b)  
{
    complex c;
    c.real = a.real*b.real - a.img*b.img;
    c.img = a.real*b.img + a.img*b.real;
    return c;
}

//複數減法的定義
complex sub(complex a, complex b)  
{
    complex c;
    c.real = a.real - b.real;
    c.img = a.img - b.img;
    return c;
}

二、第二個要解決的問題就是蝶形運算


 

　　對N個輸入行，m表示第m列蝶形。blog

　　①第1級(第1列)每一個蝶形的兩節點「距離」爲1，第2級每一個蝶形的兩節點「距離」爲2，第3級每一個蝶形的兩節點「距離」爲4，第4級每一個蝶形的兩節點「距離」爲8。由此推得，
　　第m級蝶形運算，每一個蝶形的兩節點「距離」L=2^(m-1)。
　　②對於16點的FFT，第1級有16組蝶形，每組有1個蝶形;第2級有4組蝶形，每組有2個蝶形;第3級有2組蝶形，每組有4個蝶形;第4級有1組蝶形，每組有8個蝶形。由此可推出，
　　對於N點的FFT，第m級有N/(2L)組蝶形，每組有L=2^(m-1)個蝶形。
　　③旋轉因子的肯定
　　以16點FFT爲例，第m級第k個旋轉因子爲，其中k=0～2^(m-1)-1，即第m級共有2^(m-1)個旋轉因子，根據旋轉因子的可約性，，因此第m級第k個旋轉因子爲，其中k=0～2^(m-1)-1。

　　生成一個的序列，代碼以下：ip

/*產生一個週期歐拉形式的Wn的值*/
void InitWn(complex *w,int n)   //n爲輸入的個數,w爲權值Wn
{
    int k;
    for (k = 0; k < n; k++)
    {
        w[k].real = cos(2 * PI / n * k);   //用歐拉公式計算旋轉因子
        w[k].img = -1 * sin(2 * PI / n * k);
    }
}

三、算法實現

　　咱們已經知道，N點FFT從左到右共有log2N級蝶形，每級有N/2L組，每組有L個。因此FFT的C語言編程只需用3層循環便可實現：最外層循環完成每一級的蝶形運算(整個FFT共log2N級)，中間層循環完成每一組的蝶形運算(每一級有N/2L組)，最內層循環完成單獨1個蝶形運算(每一組有L個)。
　內存

/*快速傅里葉變換*/
void fft(complex *IN_X, int n)
{
    int row = log(n) / log(2);    //row爲FFT的N個輸入對應的最大列數
    int i = 0, j = 0, k = 0, L = 0;   
    complex top, bottom, product;
    Reverse(IN_X,n);  //調用變址函數
    for (i = 0; i < row ; i++)        /*第i級蝶形運算 */
    {
        L = 1 << i;  //L等於2的i次方,表示每一級蝶形中組間的距離，即一組中蝶形個數
        for (j = 0; j < n; j += 2 * L)    /*j爲組偏移地址，每L個蝶形是一組，每級有N/2L組*/
        {
            for (k = 0; k < L; k++)        /*k爲一組中蝶形的偏移地址，一個蝶形運算 每一個group內的蝶形運算*/
            {
                product = mul(IN_X[j + k + L], Wn[n*k/2/L]);
                top = add(IN_X[j + k], product);
                bottom = sub(IN_X[j + k], product);
                IN_X[j + k] = top;
                IN_X[j + k + L] = bottom;
            }
        }
    }
}

4.所有代碼

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <malloc.h>
#define N 1000
#define PI atan(1)* 4

/*定義複數類型*/
typedef struct {
    double real;
    double img;
}complex;

complex x[N], *Wn;                 /*輸入序列,複數係數Wn*/
int size_x;                         /*size_x爲採樣的信號個數，必須爲2的倍數*/
void fft(complex *IN_X, int n);         /*對輸入IN_X，快速傅里葉變換*/
void InitWn(complex *w, int n);  /*生成一個長度爲n的Wn歐拉形式序列*/
void Reverse(complex *IN_X, int n); /*對輸入IN_X地址*/
complex add(complex, complex); /*複數加法*/
complex mul(complex, complex); /*複數乘法*/
complex sub(complex, complex); /*複數減法*/
void output();                   /*輸出快速傅里葉變換的結果*/


int main()
{
    int i;                             /*輸出結果*/
    system("cls"); //調用系統命令，CLS的功能爲清除屏幕輸出
    printf("輸出DIT方法實現的FFT結果\n");
    printf("Please input the size of x:\n");//輸入序列的大小
    scanf_s("%d", &size_x);
    printf("Please input the data in x[N]:\n");//輸入序列的實部和虛部
    for (i = 0; i < size_x; i++)
    {
        printf("請輸入第%d個序列：", i);
        scanf_s("%lf%lf", &x[i].real, &x[i].img);
    }
    printf("輸出倒序後的序列\n");
    Wn = (complex *)malloc(sizeof(complex) * size_x);  //分配變換後的值的內存空間
    InitWn(Wn , size_x);//調用變換核
    fft(x, size_x);//調用快速傅里葉變換
    printf("輸出FFT後的結果\n");
    output();//調用輸出傅里葉變換結果函數
    scanf_s("%d", &size_x);
    //free(Wn);
    return 0;
}


/*快速傅里葉變換*/
void fft(complex *IN_X, int n)
{
    int row = log(n) / log(2);    //row爲FFT的N個輸入對應的最大列數
    int i = 0, j = 0, k = 0, L = 0;   
    complex top, bottom, product;
    Reverse(IN_X,n);  //調用變址函數
    for (i = 0; i < row ; i++)        /*第i級蝶形運算 */
    {
        L = 1 << i;  //L等於2的i次方,表示每一級蝶形中組間的距離，即一組中蝶形個數
        for (j = 0; j < n; j += 2 * L)    /*j爲組偏移地址，每L個蝶形是一組，每級有N/2L組*/
        {
            for (k = 0; k < L; k++)        /*k爲一組中蝶形的偏移地址，一個蝶形運算 每一個group內的蝶形運算*/
            {
                product = mul(IN_X[j + k + L], Wn[n*k/2/L]);
                top = add(IN_X[j + k], product);
                bottom = sub(IN_X[j + k], product);
                IN_X[j + k] = top;
                IN_X[j + k + L] = bottom;
            }
        }
    }
}


/*產生一個週期歐拉形式的Wn的值*/
void InitWn(complex *w,int n)   //n爲輸入的個數,w爲權值Wn
{
    int k;
    for (k = 0; k < n; k++)
    {
        w[k].real = cos(2 * PI / n * k);   //用歐拉公式計算旋轉因子
        w[k].img = -1 * sin(2 * PI / n * k);
    }
}


/*變址計算，將x(n)碼位倒置*/
/*
碼位倒序要解決兩個問題：
    ①將t位二進制數倒序;②將倒序後的兩個存儲單元進行交換。
    若是輸入序列的天然順序號i用二進制數表示，
    例如若最大序號爲15，即用4位就可表示n3n2n1n0，
    則其倒序後j對應的二進制數就是n0n1n2n3
*/
void Reverse(complex *IN_X, int n)
{
    int row = log(n) / log(2); //row爲FFT的N個輸入對應的最大列數
    complex temp;  //臨時交換變量
    unsigned short i = 0, j = 0, k = 0;
    double t;
    for (i = 0; i < row; i++)  //從第0個序號到第N-1個序號
    {
        k = i;   //當前第i個序號
        j = 0;   //存儲倒序後的序號，先初始化爲0
        t = row;  //共移位t次
        while ((t--) > 0)    //利用按位與以及循環實現碼位顛倒
        {
            j = j << 1;
            j |= (k & 1);    //j左移一位而後加上k的最低位
            k = k >> 1;      //k右移一位，次低位變爲最低位
                             /*j爲倒序，k爲正序，
                               從k右向左的值依次移位，
                               j向左依次填入對應的k的右向位*/
        }
        if (j > i)   //若是倒序後大於原序數，就將兩個存儲單元進行交換(判斷j>i是爲了防止重複交換
        {
            temp = IN_X[i];
            IN_X[i] = IN_X[j];
            IN_X[j] = temp;
        }
    }
}



/*輸出傅里葉變換的結果*/
void output()
{
    int i;
    printf("The result are as follows：\n");
    for (i = 0; i < size_x; i++)
    {
        printf("%.4f", x[i].real);
        if (x[i].img >= 0.0001)printf("+%.4fj\n", x[i].img);
        else if (fabs(x[i].img) < 0.0001)printf("\n");
        else printf("%.4fj\n", x[i].img);
    }
}

//複數加法的定義
complex add(complex a, complex b)
{
    complex c;
    c.real = a.real + b.real;
    c.img = a.img + b.img;
    return c;
}

//複數乘法的定義
complex mul(complex a, complex b)
{
    complex c;
    c.real = a.real*b.real - a.img*b.img;
    c.img = a.real*b.img + a.img*b.real;
    return c;
}

//複數減法的定義
complex sub(complex a, complex b)
{
    complex c;
    c.real = a.real - b.real;
    c.img = a.img - b.img;
    return c;
}
相關文章
相關標籤/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公眾號
   歡迎關注本站公眾號,獲取更多信息
聯繫我們 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程