学习DIP第3天
傅里叶变换是一个非常大的话题。今天实现了下一维的DFT,兴许将完毕其它傅里叶系的算法实现和实验。
DFT公式:
-
当中e 是自然对数的底数,i是虚数单位。通常以符号表示这一变换,即
IDFT公式:
- 记为:
-
-
// // 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); }