美文网首页
FFT变换(快速傅里叶变换)

FFT变换(快速傅里叶变换)

作者: 客昂康 | 来源:发表于2021-03-06 09:59 被阅读0次

FFT.h:

//==============================================================================
//  Copyright (C) 2020 王小康. All rights reserved.
//
//  作者: 王小康
//  描述: 一维二维快速傅里叶变换(FFT)
//  日期: 2020-07-13
//
//==============================================================================

#ifndef _FFT_H_
#define _FFT_H_

typedef struct {
    float real;
    float imag;
} FFT_COMPLEX;

#ifdef __cplusplus
extern "C" {
#endif

void* FFT_init  (int bit);
void  FFT_uninit(void *fftCtx);

int FFT_1D_fft (void *fftCtx, FFT_COMPLEX *out, FFT_COMPLEX *in);
int FFT_1D_ifft(void *fftCtx, FFT_COMPLEX *out, FFT_COMPLEX *in);
int FFT_2D_fft (void *fftCtx, FFT_COMPLEX *out, FFT_COMPLEX *in);
int FFT_2D_ifft(void *fftCtx, FFT_COMPLEX *out, FFT_COMPLEX *in);

int FFT_1D_getSpectrum(void *fftCtx, float *out, FFT_COMPLEX *in);

#ifdef __cplusplus
}
#endif

#endif

FFT.c:

//==============================================================================
//  Copyright (C) 2020 王小康. All rights reserved.
//
//  作者: 王小康
//  描述: 一维二维快速傅里叶变换(FFT)
//  日期: 2020-07-13
//
//==============================================================================

#include <stdlib.h>
#include <stdint.h>
#include <math.h>
#include "FFT.h"

typedef struct {
    uint32_t  bit;
    uint32_t  len;
    float    *sin;
    float    *cos;
    uint16_t *index;
    uint8_t   buffer[0];
} FFT_CTX;

//左右翻转二进制位
static uint16_t reverseBit(uint16_t bit, uint16_t value){
    uint16_t i, ret = 0;
    for(i=0; i<bit; i++){
        ret = (ret<<1) | ((value>>i)&1);
    }
    return ret;
}

//对一行进行 FFT
static void horizontal_FFT(FFT_CTX *fftCtx, FFT_COMPLEX *out, FFT_COMPLEX *in){
    uint16_t a, b, c, d, e, f, g;
    FFT_COMPLEX tmp;
    for(a=0; a<fftCtx->len; a++){
        out[a].real = in[fftCtx->index[a]].real;
        out[a].imag = in[fftCtx->index[a]].imag;
    }
    for(a=0; a<fftCtx->bit; a++){
        b = 1 << a;
        c = 2 << a;
        d = 1 << (fftCtx->bit-1-a);
        for(e=0,f=0; e<b; e+=1,f+=d){
            for(g=e; g<fftCtx->len; g+=c){
                tmp.real = out[g+b].real * fftCtx->cos[f] + out[g+b].imag * fftCtx->sin[f];
                tmp.imag = out[g+b].imag * fftCtx->cos[f] - out[g+b].real * fftCtx->sin[f];
                out[g+b].real = out[g].real - tmp.real;
                out[g+b].imag = out[g].imag - tmp.imag;
                out[ g ].real = out[g].real + tmp.real;
                out[ g ].imag = out[g].imag + tmp.imag;
            }
        }
    }
}

//对一列进行 FFT
static void vertical_FFT(FFT_CTX *fftCtx, FFT_COMPLEX *out, FFT_COMPLEX *in){
    uint16_t a, b, c, d, e, f, g;
    FFT_COMPLEX tmp;
    for(a=0; a<fftCtx->len; a++){
        out[fftCtx->len*a].real = in[fftCtx->len*fftCtx->index[a]].real;
        out[fftCtx->len*a].imag = in[fftCtx->len*fftCtx->index[a]].imag;
    }
    for(a=0; a<fftCtx->bit; a++){
        b = 1 << a;
        c = 2 << a;
        d = 1 << (fftCtx->bit-1-a);
        for(e=0,f=0; e<b; e+=1,f+=d){
            for(g=e; g<fftCtx->len; g+=c){
                tmp.real = out[fftCtx->len*(g+b)].real * fftCtx->cos[f] + out[fftCtx->len*(g+b)].imag * fftCtx->sin[f];
                tmp.imag = out[fftCtx->len*(g+b)].imag * fftCtx->cos[f] - out[fftCtx->len*(g+b)].real * fftCtx->sin[f];
                out[fftCtx->len*(g+b)].real = out[fftCtx->len*g].real - tmp.real;
                out[fftCtx->len*(g+b)].imag = out[fftCtx->len*g].imag - tmp.imag;
                out[  fftCtx->len*g  ].real = out[fftCtx->len*g].real + tmp.real;
                out[  fftCtx->len*g  ].imag = out[fftCtx->len*g].imag + tmp.imag;
            }
        }
    }
}

//对一行进行 IFFT
static void horizontal_IFFT(FFT_CTX *fftCtx, FFT_COMPLEX *out, FFT_COMPLEX *in){
    uint16_t a, b, c, d, e, f, g;
    FFT_COMPLEX tmp;
    for(a=0; a<fftCtx->len; a++){
        out[a].real = in[fftCtx->index[a]].real;
        out[a].imag = in[fftCtx->index[a]].imag;
    }
    for(a=0; a<fftCtx->bit; a++){
        b = 1 << a;
        c = 2 << a;
        d = 1 << (fftCtx->bit-1-a);
        for(e=0,f=0; e<b; e+=1,f+=d){
            for(g=e; g<fftCtx->len; g+=c){
                tmp.real = out[g+b].real * fftCtx->cos[f] - out[g+b].imag * fftCtx->sin[f];
                tmp.imag = out[g+b].imag * fftCtx->cos[f] + out[g+b].real * fftCtx->sin[f];
                out[g+b].real = out[g].real - tmp.real;
                out[g+b].imag = out[g].imag - tmp.imag;
                out[ g ].real = out[g].real + tmp.real;
                out[ g ].imag = out[g].imag + tmp.imag;
            }
        }
    }
    for(a=0; a<fftCtx->len; a++){
        out[a].real /= fftCtx->len;
        out[a].imag /= fftCtx->len;
    }
}

//对一列进行 IFFT
static void vertical_IFFT(FFT_CTX *fftCtx, FFT_COMPLEX *out, FFT_COMPLEX *in){
    uint16_t a, b, c, d, e, f, g;
    FFT_COMPLEX tmp;
    for(a=0; a<fftCtx->len; a++){
        out[fftCtx->len*a].real = in[fftCtx->len*fftCtx->index[a]].real;
        out[fftCtx->len*a].imag = in[fftCtx->len*fftCtx->index[a]].imag;
    }
    for(a=0; a<fftCtx->bit; a++){
        b = 1 << a;
        c = 2 << a;
        d = 1 << (fftCtx->bit-1-a);
        for(e=0,f=0; e<b; e+=1,f+=d){
            for(g=e; g<fftCtx->len; g+=c){
                tmp.real = out[fftCtx->len*(g+b)].real * fftCtx->cos[f] - out[fftCtx->len*(g+b)].imag * fftCtx->sin[f];
                tmp.imag = out[fftCtx->len*(g+b)].imag * fftCtx->cos[f] + out[fftCtx->len*(g+b)].real * fftCtx->sin[f];
                out[fftCtx->len*(g+b)].real = out[fftCtx->len*g].real - tmp.real;
                out[fftCtx->len*(g+b)].imag = out[fftCtx->len*g].imag - tmp.imag;
                out[  fftCtx->len*g  ].real = out[fftCtx->len*g].real + tmp.real;
                out[  fftCtx->len*g  ].imag = out[fftCtx->len*g].imag + tmp.imag;
            }
        }
    }
    for(a=0; a<fftCtx->len; a++){
        out[fftCtx->len*a].real /= fftCtx->len;
        out[fftCtx->len*a].imag /= fftCtx->len;
    }
}

////////////////////////////////////////////////////////////////////////////////////////////////////

//初始化，生成变址表、sin表和cos表。
void* FFT_init(int bit){
    if((bit < 3) || (bit > 15)) return NULL;
    uint32_t len = 1 << bit;
    FFT_CTX *ctx = malloc(sizeof(FFT_CTX) + (sizeof(float)+sizeof(uint16_t))*len);
    if(ctx == NULL) return NULL;
    ctx->bit   = bit;
    ctx->len   = len;
    ctx->sin   = (float*)(ctx->buffer) + 0;
    ctx->cos   = (float*)(ctx->buffer) + len/2;
    ctx->index = (uint16_t*)((float*)(ctx->buffer) + len);
    uint32_t i;
    for(i=0; i<len; i++){
        ctx->index[i] = reverseBit(bit&0x00ff, i);
    }
    for(i=0; i<len/2; i++){
        ctx->sin[i] = sin((3.14159265359*2/len)*i);
        ctx->cos[i] = cos((3.14159265359*2/len)*i);
    }
    return ctx;
}

//去初始化
void FFT_uninit(void *fftCtx){
    free(fftCtx);
}

//一维 FFT
int FFT_1D_fft(void *fftCtx, FFT_COMPLEX *out, FFT_COMPLEX *in){
    horizontal_FFT(fftCtx, out, in);
    return 0;
}

//一维 IFFT
int FFT_1D_ifft(void *fftCtx, FFT_COMPLEX *out, FFT_COMPLEX *in){
    horizontal_IFFT(fftCtx, out, in);
    return 0;
}

//二维 FFT。先一行行做 FFT，再一列列做 FFT。
int FFT_2D_fft(void *fftCtx, FFT_COMPLEX *out, FFT_COMPLEX *in){
    uint32_t i;
    FFT_CTX *ctx = fftCtx;
    FFT_COMPLEX *tmp = malloc(ctx->len*ctx->len*sizeof(FFT_COMPLEX));
    if(tmp == NULL) return -1;
    for(i=0; i<ctx->len; i++){
        horizontal_FFT(ctx, tmp+(i*ctx->len), in+(i*ctx->len));
    }
    for(i=0; i<ctx->len; i++){
        vertical_FFT(ctx, out+i, tmp+i);
    }
    free(tmp);
    return 0;
}

//二维 IFFT。先一列列做 IFFT，再一行行做 IFFT。
int FFT_2D_ifft(void *fftCtx, FFT_COMPLEX *out, FFT_COMPLEX *in){
    uint32_t i;
    FFT_CTX *ctx = fftCtx;
    FFT_COMPLEX *tmp = malloc(ctx->len*ctx->len*sizeof(FFT_COMPLEX));
    if(tmp == NULL) return -1;
    for(i=0; i<ctx->len; i++){
        vertical_IFFT(ctx, tmp+i, in+i);
    }
    for(i=0; i<ctx->len; i++){
        horizontal_IFFT(ctx, out+(i*ctx->len), tmp+(i*ctx->len));
    }
    free(tmp);
    return 0;
}

//获取一维FFT后的频谱图
int FFT_1D_getSpectrum(void *fftCtx, float *out, FFT_COMPLEX *in){
    uint32_t i;
    FFT_CTX *ctx = fftCtx;
    out[0] = sqrt(in[0].real*in[0].real+in[0].imag*in[0].imag) / ctx->len;
    for(i=1; i<=ctx->len/2; i++){
        out[i] = 2 * sqrt(in[i].real*in[i].real + in[i].imag*in[i].imag) / ctx->len;
    }
    return 0;
}

测试程序test.c:

//==============================================================================
//  Copyright (C) 2020 王小康. All rights reserved.
//
//  作者: 王小康
//  描述: 一维二维快速傅里叶变换(FFT)测试程序
//  日期: 2020-07-13
//
//==============================================================================

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include "FFT.h"

#define FFT_2PI     (3.14159265359*2)
#define FFT_1D_BIT  (8)
#define FFT_1D_LEN  (1<<FFT_1D_BIT)
#define FFT_2D_BIT  (4)
#define FFT_2D_LEN  (1<<FFT_2D_BIT)

static void showWave(FFT_COMPLEX *wave, int num, float ratio){
    int i, j, A;
    for(i=0; i<num; i++){
        printf("%13f %13f ", wave[i].real, wave[i].imag);
        A = (int)(wave[i].real*ratio+(wave[i].real>=0.0?0.5:-0.5)) + 80;
        for(j=0; j<=160; j++){
            if     (j == A ) putchar('0');
            else if(j == 80) putchar('|');
            else             putchar('-');
        }
        putchar('\n');
    }
}

static void showSpectrum(float *wave, int num){
    int i, j, A;
    for(i=0; i<num; i++){
        printf("%11dHz %13f ", i, wave[i]);
        A = (int)(wave[i]+0.5);
        for(j=0; j<=160; j++){
            if(j < A ) putchar('|');
            else       putchar('-');
        }
        putchar('\n');
    }
}

static void test1D(void){
    float buffer0[1+FFT_1D_LEN/2];
    FFT_COMPLEX buffer1[FFT_1D_LEN];
    FFT_COMPLEX buffer2[FFT_1D_LEN];
    FFT_COMPLEX buffer3[FFT_1D_LEN];
    void *fftCtx = FFT_init(FFT_1D_BIT);
    
    //通过多个正弦波构造复合波，FFT_1D_LEN 个采样点视为一个周期。
    int i;
    for(i=0; i<FFT_1D_LEN; i++){
        buffer1[i].real = 10;                                           //直流分量 10
        buffer1[i].real += sin(i*(FFT_2PI/FFT_1D_LEN) *  1 + 10) * 30;  //频率1Hz， 振幅30，初相位10。
        buffer1[i].real += sin(i*(FFT_2PI/FFT_1D_LEN) *  2 + 30) * 40;  //频率2Hz， 振幅40，初相位30。
        buffer1[i].real += sin(i*(FFT_2PI/FFT_1D_LEN) *  3 + 20) * 20;  //频率3Hz， 振幅20，初相位20。
        buffer1[i].real += sin(i*(FFT_2PI/FFT_1D_LEN) *  6 + 30) * 40;  //频率6Hz， 振幅40，初相位30。
        buffer1[i].real += sin(i*(FFT_2PI/FFT_1D_LEN) *  9 + 15) * 16;  //频率9Hz， 振幅16，初相位15。
        buffer1[i].real += sin(i*(FFT_2PI/FFT_1D_LEN) * 10 + 13) * 15;  //频率10Hz，振幅15，初相位13。
        buffer1[i].real += sin(i*(FFT_2PI/FFT_1D_LEN) * 12 + 11) * 11;  //以此类推
        buffer1[i].real += sin(i*(FFT_2PI/FFT_1D_LEN) * 15 + 25) * 13;  //以此类推
        buffer1[i].real += sin(i*(FFT_2PI/FFT_1D_LEN) * 16 + 15) * 20;  //以此类推
        buffer1[i].real += sin(i*(FFT_2PI/FFT_1D_LEN) * 19 + 19) * 10;  //以此类推
        buffer1[i].real += sin(i*(FFT_2PI/FFT_1D_LEN) * 20 + 14) * 16;  //以此类推
        buffer1[i].real += sin(i*(FFT_2PI/FFT_1D_LEN) * 25 + 21) * 10;  //以此类推
        buffer1[i].imag = 0;
    }
    printf("原始波:\n");
    showWave(buffer1, FFT_1D_LEN, 0.35);
    putchar('\n');
    
    //快速傅里叶变换 获取频谱
    FFT_1D_fft(fftCtx, buffer2, buffer1);
    FFT_1D_getSpectrum(fftCtx, buffer0, buffer2);
    printf("进行FFT后获取的频谱:\n");
    showSpectrum(buffer0, 1+FFT_1D_LEN/2);
    putchar('\n');
    
    //快速傅里叶逆变换
    FFT_1D_ifft(fftCtx, buffer3, buffer2);
    printf("进行IFFT之后:\n");
    showWave(buffer3, FFT_1D_LEN, 0.35);
    putchar('\n');
    
    FFT_uninit(fftCtx);
}

void show2D(FFT_COMPLEX data[][FFT_2D_LEN]){
    int i, j;
    for(i=0; i<FFT_2D_LEN; i++){
        for(j=0; j<FFT_2D_LEN; j++) printf("%13f", data[i][j].real); putchar('\n');
        for(j=0; j<FFT_2D_LEN; j++) printf("%13f", data[i][j].imag); putchar('\n');
        putchar('\n');
    }
    putchar('\n');
}

static void test2D(void){
    FFT_COMPLEX buffer1[FFT_2D_LEN][FFT_2D_LEN];
    FFT_COMPLEX buffer2[FFT_2D_LEN][FFT_2D_LEN];
    FFT_COMPLEX buffer3[FFT_2D_LEN][FFT_2D_LEN];
    void *fftCtx = FFT_init(FFT_2D_BIT);
    
    //构造数据
    int i, j;
    for(i=0; i<FFT_2D_LEN; i++){
        for(j=0; j<FFT_2D_LEN; j++){
            buffer1[i][j].real = i;
            buffer1[i][j].real += sin(j*(FFT_2PI/FFT_2D_LEN) *  1 + 10 + i) * 30;
            buffer1[i][j].real += sin(j*(FFT_2PI/FFT_2D_LEN) *  2 + 15 + i) * 40;
            buffer1[i][j].real += sin(j*(FFT_2PI/FFT_2D_LEN) *  3 + 20 + i) * 20;
            buffer1[i][j].real += sin(j*(FFT_2PI/FFT_2D_LEN) *  6 + 30 + i) * 40;
            buffer1[i][j].imag = 0;
        }
    }
    printf("原始矩阵:\n");
    show2D(buffer1);
    
    //2维 FFT
    FFT_2D_fft(fftCtx, (void*)buffer2, (void*)buffer1);
    printf("进行二维FFT之后的矩阵:\n");
    show2D(buffer2);
    
    //2维 IFFT
    FFT_2D_ifft(fftCtx, (void*)buffer3, (void*)buffer2);
    printf("进行二维IFFT之后的矩阵:\n");
    show2D(buffer3);
    
    FFT_uninit(fftCtx);
}

static void help(char *name){
    printf("Usage: %s 1|2\n", name);
    printf("       1: 1D FFT test.\n");
    printf("       2: 2D FFT test.\n");
}

int main(int argc, char* argv[]){
    if(argc < 2){
        help(argv[0]);
        return 0;
    }
    if(strcmp(argv[1], "1") == 0){
        test1D();
    }
    else if(strcmp(argv[1], "2") == 0){
        test2D();
    }
    else {
        help(argv[0]);
        return 0;
    }
    return 0;
}

测试结果


原始波形(局部)
FFT变换后获取的频谱(局部)
FFT逆变换后的波形(局部)

相关文章

网友评论

      本文标题：FFT变换(快速傅里叶变换)

      本文链接：https://www.haomeiwen.com/subject/puujqltx.html

      延伸阅读

      深度阅读

      • 您也可以注册成为美文阅读网的作者，发表您的原创作品、分享您的心情！

      栏目导航

      热点阅读

      关于我们|服务条款|联系我们|FFT变换(快速傅里叶变换)|投稿指南|网站地图|RSS订阅|排版工具|手机版

      提供经典美文摘抄,优美散文欣赏,现代诗歌精选,短篇小说,心情随笔,表白情书范文,故事会在线阅读欣赏

      Copyright © 2014-2023 Haomeiwen.com All Rights Reserved. 好美文阅读网 版权所有

      备案信息：桂公网安备 45052102000051号 · 桂ICP备13007215号-3

      本站所收录作品、热点评论等信息部分来源互联网，目的只是为了系统归纳学习和传递资讯

      所有作品版权归原创作者所有，与本站立场无关，如不慎侵犯了你的权益，请联系我们告知，我们将做删除处理！