zoukankan      html  css  js  c++  java
  • 快速卷积

    一、功能

    用快速傅里叶变换计算两个有限长序列的快速卷积。

    二、方法简介

    设序列(x(n))的长度为(M),序列(y(n))的长度为(N),序列(x(n))(y(n))的线性卷积定义为

    [z(n)=sum_{i=0}^{M-1}x(i)y(n-i) , n=0,1,...,M+N-2 ]

    用快速傅里叶变换计算线性卷积的算法如下

    1、选择(L)满足下述条件

    [left{egin{matrix}egin{align*}L &geqslant M + N - 1\ L &= 2^{gamma }, gamma is a positive integerend{align*}end{matrix} ight. ]

    2、将序列(x(n))(y(n))按如下方式补零,形成长为(L = 2^{gamma })的序列

    [egin{matrix}x(n)=left{egin{matrix}egin{align*}x(n) &, n=0,1,...,M-1 \ 0 &, n=M,M+1,...,L-1end{align*}end{matrix} ight.\ end{matrix} ]

    [egin{matrix}y(n)=left{egin{matrix}egin{align*}y(n) &, n=0,1,...,N-1 \ 0 &, n=N,N+1,...,L-1end{align*}end{matrix} ight.\ end{matrix} ]

    3、用FFT算法分别计算(x(n))(y(n))的离散傅里叶变换(X(k))(Y(k))

    [egin{matrix}X(k)=sum_{n=0}^{L-1}x(n)e^{-j2pi nk/L}\ Y(k)=sum_{n=0}^{L-1}y(n)e^{-j2pi nk/L}end{matrix} ]

    4、计算(X(k))(Y(k))的乘积

    [Z(k)=X(k)Y(K) ]

    5、用FFT算法计算(Z(k))的离散傅里叶反变换,得到卷积(z(n))

    [z(n)=frac{1}{L}sum_{k=0}^{L-1}Z(k)e^{j2pi nk/L}, n=0,1,...,L-1 ]

    序列(z(n))的前(M+N-1)点的值就是序列(x(n))(y(n))的线性卷积。

    三、使用说明

    快速卷积的C语言实现方式如下

    /************************************
    	x		----双精度一维数组,长度为len。开始时存放实序列x(i),最后存放线性卷积的值。
    	y		----双精度一维数组,长度为n。开始时存放实序列y(i)。
    	m		----数据长度,序列x(i)的长度。
    	n		----数据长度,序列y(i)的长度。
    	len		----线性卷积长度,len≥m+n-1,且必须是2的整数次幂,即len=2^gamma。
    ************************************/
    #include "rfft.c"
    #include "irfft.c"
    void convol(double *x, double *y, int m, int n, int len)
    {
    	int i, len2;
    	double t;
    	for(i = m; i < len; i++)
    		x[i] = 0.0;
    	for(i = n; i < len; i++)
    		y[i] = 0.0;
    	rfft(x, len);
    	irfft(y, len);
    	len2 = len / 2;
    	x[0] = x[0] * y[0];
    	x[len2] = x[len2] * y[len2];
    	for( i = 1; i < len2; i++){
    		t = x[i] * y[i] - x[len - i] * y[len - i];
    		x[len - i] = x[i] * y[len - i] + x[len - i] * y[i];
    		x[i] = t;
    	}
    	irfft(x, len);
    }
    

    其中rfft.c文件请参考实序列快速傅里叶变换(一)

    irfft.c在rfft.c的基础上添加系数长度的倒数。

  • 相关阅读:
    自学Linux命令的四种方法
    POJ 1170 Shopping Offers -- 动态规划(虐心的六重循环啊!!!)
    九度OJ 1447 最短路 1008 最短路径问题
    九度OJ 1024 畅通工程 -- 并查集、贪心算法(最小生成树)
    PHPActiveRecord 学习三
    PHPUnit 组织测试
    PHPActiveRecord validates
    PHPActiveRecord 学习二
    PHPActiveRecord 学习一
    PHP ActiveRecord demo栗子中 关于类名 的问题
  • 原文地址:https://www.cnblogs.com/liam-ji/p/11979964.html
Copyright © 2011-2022 走看看