zoukankan      html  css  js  c++  java
  • 一维DFT

    学习DIP第3天

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

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

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


           IDFT公式:

                                     xleft[n
ight]={1 over N}sum_{k=0}^{N-1} e^{ ifrac{2pi}{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:
    ");
        srand((int)time(0));
        for(int i=0;i<SIZE;i++){
            data[i]=rand()%VALUE_MAX;
            printf("%lf
    ",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
    ",real,imagin);
            else
                printf("%lf%lfj
    ",real,imagin);
        }
        end=clock();
        printf("DFT use time :%lf for Datasize of:%d
    ",(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
    ",real,imagin);
            else
                printf("%lf%lfj
    ",real,imagin);
        }
        end=clock();
        printf("IDFT use time :%lfs for Datasize of:%d
    ",(double)(end-start)/CLOCKS_PER_SEC,size);
        
    
    
    }
    int main(int argc, const char * argv[]) {
        double input[SIZE];
        Complex dst[SIZE];
        setInput(input,SIZE);
        printf("
    
    ");
        DFT(input, dst, SIZE);
        printf("
    
    ");
        IDFT(dst, NULL, SIZE);
        
    }
    



  • 相关阅读:
    sql语句查询数据库中含有某字符串的表名
    PHP复制文件夹及文件夹内的文件
    Vue.js绑定内联样式
    Vue模板语法V-bind
    Vue实例
    Vue.js几个简单用法
    Git Pull Failed: cannot lock ref 'refs/remotes/origin/xxxxxxxx': unable to resolve ref
    SSM 框架详细整合教程(IDEA版)(Spring+SpringMVC+MyBatis)
    IntelliJ IDEA(2017/2018)安装图解与破解教程
    Hadoop集群单机伪分布搭建
  • 原文地址:https://www.cnblogs.com/mengfanrong/p/5239836.html
Copyright © 2011-2022 走看看