基於fftw的二維FFT實現

作者使用的fftw版本為3.3.5

fftw官網下載地址http://www.fftw.org/install/windows.html

#include "fftw3.h"
#include <stdio.h>
#include <stdlib.h>
#include <opencv2/opencv.hpp>
#include <iostream>

using namespace std;
using namespace cv;

typedef struct _COMPLEX
{
    float real;
    float img;
}COMPLEX;

float **m;
COMPLEX **fm;
float **ffm;

//申請動態內存
float** allocfloat(float **mem, int h, int w)
{
    mem = (float **)malloc(sizeof(float *) * h);//分配指針數組
    mem[0] = (float *)malloc(sizeof(float) * h * w);//一次性分配所有空間
    for (int i = 1; i<h; i++)//注意從1開始
    {
        mem[i] = mem[i - 1] + w;
    }
    return(mem);
}
//申請動態內存
COMPLEX** alloccomplex(COMPLEX **mem, int h, int w)
{
    mem = (COMPLEX **)malloc(sizeof(COMPLEX *) * h);//分配指針數組
    mem[0] = (COMPLEX *)malloc(sizeof(COMPLEX) * h * w);//一次性分配所有空間
    for (int i = 1; i<h; i++)//注意從1開始
    {
        mem[i] = mem[i - 1] + w;
    }
    return(mem);
}

void FFT2(float **input, COMPLEX **output, int height, int width)
{
    fftwf_plan planR;
    fftwf_complex *inR, *outR;
    inR = (fftwf_complex*)fftw_malloc(sizeof(fftwf_complex) * height * width);
    outR = (fftwf_complex*)fftw_malloc(sizeof(fftwf_complex) * height * width);
    for (int i = 0; i < height; i++)
        for (int j = 0; j < width; j++)
        {
            inR[i * width + j][0] = input[i][j];
            inR[i * width + j][1] = 0;
        }
    planR = fftwf_plan_dft_2d(height, width, inR, outR, FFTW_FORWARD, FFTW_ESTIMATE);
    fftwf_execute(planR);
    for (int i = 0; i < height; i++)
        for (int j = 0; j < width; j++)
        {
            output[i][j]._COMPLEX::real = outR[i * width + j][0];
            output[i][j]._COMPLEX::img = outR[i * width + j][1];
        }
    fftwf_destroy_plan(planR);
    fftwf_free(inR);
    fftwf_free(outR);
}

void IFFT2(COMPLEX **input, float **output, int height, int width)
{
    fftwf_plan planR;
    fftwf_complex *inR, *outR;
    inR = (fftwf_complex*)fftw_malloc(sizeof(fftwf_complex) * height * width);
    outR = (fftwf_complex*)fftw_malloc(sizeof(fftwf_complex) * height * width);
    for (int i = 0; i < height; i++)
        for (int j = 0; j < width; j++)
        {
            inR[i * width + j][0] = input[i][j]._COMPLEX::real;
            inR[i * width + j][1] = input[i][j]._COMPLEX::img;
        }
    planR = fftwf_plan_dft_2d(height, width, inR, outR, FFTW_BACKWARD, FFTW_ESTIMATE);
    fftwf_execute(planR);
    for (int i = 0; i < height; i++)
        for (int j = 0; j < width; j++)
        {
            output[i][j] = outR[i * width + j][0] / (height*width);
        }
    fftwf_destroy_plan(planR);
    fftwf_free(inR);
    fftwf_free(outR);
}

void main()
{int width = 6;
    int height = 10;
    m = allocfloat(m, height, width);
    ffm = allocfloat(ffm, height, width);
    fm = alloccomplex(fm, height, width);
    for (int i = 0; i < height; i++)
    {
        for (int j = 0; j < width; j++)
        {
            m[i][j] = i*width + j;

        }
    }
    for (int i = 0; i < height; i++)
    {
        for (int j = 0; j < width; j++)
        {
            printf("%.2f ", m[i][j]);
        }
        printf("\n");
    }
    FFT2(m, fm, height, width);
    IFFT2(fm, ffm, height, width);
    for (int i = 0; i < height; i++)
    {
        for (int j = 0; j < width; j++)
        {
            printf("%.2f ", ffm[i][j]);
        }
        printf("\n");
    }
    while (1);
}