一维DFT

学习DIP第3天

       傅里叶变换是一个非常大的话题。今天实现了下一维的DFT,兴许将完毕其它傅里叶系的算法实现和实验。

       DFT公式:
                       \hat{x}[k]=\sum_{n=0}^{N-1} e^{-i\frac{2\pi}{N}nk}x[n] \qquad k = 0,1,\ldots,N-1.

    当中e 是自然对数的底数i是虚数单位。通常以符号\mathcal{F}表示这一变换,即
                 \hat{x}=\mathcal{F}x


       IDFT公式:

                                  x\left[n\right]={1 \over N}\sum_{k=0}^{N-1} e^{ i\frac{2\pi}{N}nk}\hat{x}[k] \qquad n = 0,1,\ldots,N-1.
   记为:
                             x=\mathcal{F}^{-1}\hat{x}
   c语言代码:
//
//  main.c
//  Fourer1D
//
//  Created by Tony on 14/11/16.
//  Copyright (c) 2014年 Tony. All rights reserved.
//

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#define SIZE 1000
#define VALUE_MAX 2000
struct Complex_{
    double real;
    double imagin;
};
typedef struct Complex_ Complex;

void setInput(double * data,int n){
    printf("Setinput signal:\n");
    srand((int)time(0));
    for(int i=0;i<SIZE;i++){
        data[i]=rand()%VALUE_MAX;
        printf("%lf\n",data[i]);
    }

}
void DFT(double * src,Complex * dst,int size){
    clock_t start,end;
    start=clock();
    
    for(int m=0;m<size;m++){
        double real=0.0;
        double imagin=0.0;
        for(int n=0;n<size;n++){
            double x=M_PI*2*m*n;
            real+=src[n]*cos(x/size);
            imagin+=src[n]*(-sin(x/size));
        
        }
        dst[m].imagin=imagin;
        dst[m].real=real;
        if(imagin>=0.0)
            printf("%lf+%lfj\n",real,imagin);
        else
            printf("%lf%lfj\n",real,imagin);
    }
    end=clock();
    printf("DFT use time :%lf for Datasize of:%d\n",(double)(end-start)/CLOCKS_PER_SEC,size);

}

void IDFT(Complex *src,Complex *dst,int size){
    //Complex temp[SIZE];
    clock_t start,end;
    start=clock();
    for(int m=0;m<size;m++){
        double real=0.0;
        double imagin=0.0;
        for(int n=0;n<size;n++){
            double x=M_PI*2*m*n/size;
            real+=src[n].real*cos(x)-src[n].imagin*sin(x);
            imagin+=src[n].real*sin(x)+src[n].imagin*cos(x);
            
        }
        real/=SIZE;
        imagin/=SIZE;
        if(dst!=NULL){
            dst[m].real=real;
            dst[m].imagin=imagin;
        }
        if(imagin>=0.0)
            printf("%lf+%lfj\n",real,imagin);
        else
            printf("%lf%lfj\n",real,imagin);
    }
    end=clock();
    printf("IDFT use time :%lfs for Datasize of:%d\n",(double)(end-start)/CLOCKS_PER_SEC,size);
    


}
int main(int argc, const char * argv[]) {
    double input[SIZE];
    Complex dst[SIZE];
    setInput(input,SIZE);
    printf("\n\n");
    DFT(input, dst, SIZE);
    printf("\n\n");
    IDFT(dst, NULL, SIZE);
    
}



评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值