请教:C或C++中卷积的快速算法

请教:C或C++中卷积的快速算法,第1张

卷积在工程和数学上都有很多应用:

统计学中,加权的滑动平均是一种卷积。概率论中,两个统计独立变量X与Y的和的概率密度函数是X与Y的概率密度函数的卷积。声学中,回声可以用源声与一个反映各种反射效应的函数的卷积表示。电子工程与信号处理中,任一个线性系统的输出都可以通过将输入信号与系统函数(系统的冲激响应)做卷积获得。物理学中,任何一个线性系统(符合叠加原理)都存在卷积。

介绍一个实际的概率学应用例子。假设需求到位时间的到达率为poisson(λ)分布,需求的大小的分布函数为D(.),则单位时间的需求量的分布函数为 F(x):

其中 D(k)(x)为k阶卷积。

卷积是一种线性运算,图像处理中常见的mask运算都是卷积,广泛应用于图谈没像滤波。castlman的书对卷积讲得很详细。

高斯变换就是用高斯函数对图像进行卷积。高斯算子可以直接从离散高斯函数得到含唯纳:

for(i=0i<Ni++)

{

for(j=0j<Nj++)

{

g[i*N+j]=exp(-((i-(N-1)/2)^2+(j-(N-1)/2)^2))/(2*delta^2))

sum += g[i*N+j]

}

}

再除以 sum 得到归一化算子

N是滤波器的大小,delta自选

首先,在提到卷积之前,必须提到卷积出现的背景。卷积是在信号与线性系统的基础上或背景中出现的,脱离这个背景单独谈卷积是没有任何意义的,除了那个所谓褶反公式上的数学意义和积分(或求和,离散情况下)。

信号与线性系统,讨论的就是信号经过一个线性系统以后发生的变化(就是输入 输出 和所经过的所谓系统,这三者之间的数学关系)。所谓线性系统的含义,就是,这个所谓的系统,带来的输出信号与输入信号的数学关系式之间是线性的运算关系。

因此,实际上,都是要根据我们需要待处理的信号形式,来设计所谓的系统传递函数,那么这个系统的传递函数和输入信号,在数学上的形式就是所谓的卷积关系。

卷积关系最重要的一种情况,就是在信号与线性系统或数字信号处理中的卷积定理。利用该定理,可以将时间域或空间域中的卷积运算等价为频率域的相乘运山大算,从而利用FFT等快速算法,实现有效的计算,节省运算代价。

C++语言代码: void convolution(float *input1, float *input2, float *output, int mm, int nn){ float *xx = new float[mm+nn-1]// do convolution for (int i = 0i <mm+nn-1i++){ xx[i] = 0.0 for (int j = 0j <mmj++) { if (i-j >0 &&i-j <nn)xx[i] += input1[j] * input2[i-j]} } // set value to the output array for (int i = 0i <mmi++) output[i] = xx[i + (nn-1) / 2]delete[] xx}

一、实验目的

1、学会FFT算法程序(或函数)的使用方法

2、了解序列的线性卷积和圆周卷积之间的关系

3、验证有限长FFT算法实现线性相关运算的快速计算方法

4、解FFT的点数对快速卷积与快速相关运算结果所产生的影响

5、了学会利用FFT算法进行有限长序列的线性卷积的快速计算

6、掌握基-2快速傅立叶变换(Fast FourierTransform,FFT)的算法原理及其程序实现方法.

二、实验原理

1、有限长序列的线性卷积和圆周卷积线性卷积和圆周卷对于有限长序列,存在两种形式的销则卷积,即积。设x(n)是长度光JM的有限长序列,y(n)是长度为N的有限长序列,则二者的线性卷积可表示为:

M-1

线性卷积的结果序列f(n)是一个长度为L=N +M -1的有限有限长序列,且长序列。将x(n)及y(n)均补零增长为L点的二者的L点的圆周卷积可表示为:

圆周卷积的结果序列f(n)是一个长度为L的有限长序列,由圆周卷积的点数所决定。有限长序列的线性卷积和圆周卷积之间的关系可用公式表示如下:

即:圆周卷积是亏指棚线性卷积以圆周卷积的点数几为周期进行周期延拓后所取的主值序列。因而,在圆周卷积的点数大于或等于线性卷积的长度时、圆周卷积结果和线性卷积结果相等,这也是快速卷积算法的理论基础之一。

2、离散傅里叶变换的卷积性质

离散傅里叶变换的卷积性质也是快速卷积算法的另一理论基础。若f(n)是有限长序列x(n)和有限长序列y(n)的L点圆周卷积,即公式(5-2),则 f(n)的L点离散傅里叶变换为: F_c (k)=X(k)Y(k)

3、Matlab中FFT与IFFT的实现

离散傅立叶变换(Discrete Fourier Transform, DFT)实现了频域的离散化,方便了计算机处理,在数字信号处理中有着非常重要的作用。但直接计算DFT的运算量与变换长度N的平方成正比,计算量太大。而快速逗烂傅立叶变换FFT则是快速计算DFT的有效算法,大大提高了DFT的运算效率,在信号频谱的分析、滤波器频率响应的计算,以及线性卷积的快速计算等方面起着非常重要的作用。FFT 采用分组计算的方式进行DFT的快速计算,具体算法原理参看教材,在附录B中也给出了常用的基-2时间抽取FFT算法和分裂基FFT 算法的C语言程序。相应的,IFFT 则为离散傅里叶反变换,即 IDFT 的快速计算方法。在Matlab中,提供了f(t)和 ifi(t)两个函数来分别实现快速傅立叶变换的正变换和反变换。Ft(t)和if(t)两个函数是用机器语言而不是Matlab 指令写成的,执行速度很快。除了输入、输出参数的含义不同之外,这两个函数的调用方法完全相同,因此以ff()函数为例说明二者的使用方法。ff()函数常用调用格式有两种:

(1)Xk=fft(xn)

其中,xn为输入时域序列x(n),返回结果xk为x(n)的离散傅里叶变换X(k)。当xn是矩阵时(对应于多通道信号),计算xn中每一列信号的离散傅里叶变换。当xn的长度是2的整数幂,采

用基2快速算法计算,否则采用较慢的混合基算法进行计算。

(2)Xk=fft(xn,NFFT)

这种调用格式相比较于上一种调用格式,多了一个输入参数NFFT,用于指定FFT的点数。当NFFT的值是2的整数幂,采用基-2快速算法计算,否则采用较慢的混合基算法进行计算。当xn的长度大于NFFT时,对 xn进行自动截断当xn的长度小于NFFT时,在xn后自动进行补零。

4、快速卷积基本算法的原理

利用FFT进行有限长序列的线性卷积的快速计算即为快速卷积算法。按照上述相关原理,利用FFT计算线性卷积,即计算公式中的f(n)的算法步骤可用下图来表示:

其中,FFT运算的点数L应满足L≥N+ M-1。为了采用基-2的算法,常常需要L还应满足L=2M

5、快速相关算法原理

利用 FFT 讲行有限长序列相关运算的快速实现则称为快速相关算法。若(n)是长度头M的有限长序列,y(n)是长度为N的有限长序列,则一者的线性百相关序列R(m)与二者的线性卷积运算之间的关系可表示

由于线性卷积可以用FFT讲行快速计算,则按公式,相关运算也可以利用fft进行快速计算,在此不再赘述其原理。需要时,根据傅立叶变换的反折和共轭特性可以减少一次FFT的使用提高计算效率。

三、实验步骤、数据记录及处理

本实验利用Matlab中提供的fft()和 ifft()函数进行快速卷积算法和快速相关算法性卷积和圆周卷积之间的关系进行验证。具体的实验内容和实验的实现,并对有限长序列线

步骤如下所示:

其中:

1、用 Matlab生成两个有限长序列x(n) y(n)

(1)基于fft()和 ifft()函数,编程利用4点快速卷积算法计算有限长序列x(n)与y(n)的卷积,结果令为c1(n)。

(2)基于fft()和 ifft()函数,编程利用速卷积算法计算有限长序列x(n)与y(n)的卷积,结果令为c2(n)

(3)调用conv()函数计算有限长序列x(n)与y(n)的卷积,结果令为c3(n)。分别绘制序列x(n)、y(n)、c1(n)、c2(n)和c3(n)的图形。对结果进行分析,并通过实验结果验证有限长序列线生卷积和圆周卷积之间的关系。

2、设两个有限长序列x(n)和h(n)分别为:

(1)x(n)=(sin 0.4n)·R,(n+1)

(2)h(n)=(0.9)"R,o(n+1)

按图所示的快速卷积算法原理编写完整的快速卷积算法程序计算y(n)=x(n)*h(n),绘制序列图形,并通过 conv()函数对结果进行验证。

3、将实验1中实验内容4所给定的信号利用快速相关算法进行自相关序列的计算,并将结果与实验1中的结果进行对比验证。

实验程序:

clcclear allclose all

nx=0:1:2x=[1,1,1]%确定x(n)与y(n)的自变量取值范围

ny=0:1:2y=[2,2,2]%生成x(n)与h(n)

L=pow2(nextpow2(length(x)+length(y)-1))%确定FFT快速卷积的点数

xk=fft(x,L)%计算x的L点FFT,结果为x(k)

yk=fft(y,L)%计算y的L点FFT,结果为y(k)

YK=xk.*yk%计算y(k)

c1n=ifft(YK,4)%四点

c2n=ifft(YK,L)%八点

c3n=conv(x,y)%线性卷积

n1=0:1:3%确定n1的取值范围

n2=0:1:7%确定n2的取值范围

n3=0:1:4%确定n3的取值范围

%绘制x(n)与y(n)波形

figure('name','x(n)与y(n)序列')

subplot(2,1,1)stem(nx,x)grid onxlabel('n')ylabel('x(n)')title('序列x(n)')

subplot(2,1,2)stem(ny,y)grid onxlabel('n')ylabel('y(n)')title('序列y(n)')

%绘制x(n)与y(n)序列卷积波形

figure('name','x(n)与y(n)序列卷积')

subplot(3,1,1)stem(n1,c1n)grid onxlabel('n')ylabel('c1(n)')title('x(n)与y(n)的4点fft')

subplot(3,1,2)stem(n2,c2n)grid onxlabel('n')ylabel('c2(n)')title('x(n)与y(n)的8点fft')

subplot(3,1,3)stem(n3,c3n)grid onxlabel('n')ylabel('c3(n)')title('x(n)与y(n)的线性卷积')

登录后复制

clcclear allclose all

nxn=1:15nhn=1:20%确定x(n)与y(n)的自变量取值范围

xn=sin(0.4*nxn)hn=0.9.^nhn%生成x(n)与h(n)

L=pow2(nextpow2(length(xn)+length(hn)-1))%确定FFT快速卷积的点数

Xk=fft(xn,L)%计算x的L点FFT,结果为x(k)

Hk=fft(hn,L)%计算y的L点FFT,结果为y(k)

Yk=Xk.*Hk%计算YK

y1n=ifft(Yk,L)%对YK调用IFFT,求得y1(n)

y2n=conv(xn,hn)%计算y2(n)的卷积

n1=(nxn(1)+nhn(1)):1:(L+nxn(1)+nhn(1)-1)%确定n1的自变量取值范围

n2=(nxn(1)+nhn(1)):1:(nxn(15)+nhn(20))%确定n2的自变量取值范围

%绘制有限长序列想x(n)与h(n)

figure('name','x(n)与h(n)序列')

subplot(2,2,1)stem(nxn,xn)grid onxlabel('n')ylabel('x(n)')title('x(n)=(sin0.4n)*R15(n+1)')

subplot(2,2,2)stem(nhn,hn)grid onxlabel('n')ylabel('h(n)')title('h(n)=((0.9)^n)R20(n+1)')

subplot(2,2,3)stem(n1,y1n)grid onxlabel('n')ylabel('y(n)')title('快速卷积法计算y(n)')

subplot(2,2,4)stem(n2,y2n)grid onxlabel('n')ylabel('y(n)')title('conv()函数计算y(n)')

登录后复制

clcclear allclose all

n1=1:100n2=1:100%确定n1,n2取值范围

x=[101,82,66,35,31,7,20,92,154,125,85,68,38,23,10,24,83,132,131,118,90,67,60,47,41,21,16,6,4,7,14,34,45,43,48,42,28,10,8,2,0,1,5,12,14,35,46,41,30,24,16,7,4,2,8,17,36,50,62,67,71,48,28,8,13,57,122,138,103,86,63,37,24,11,15,40,62,98,124,96,66,64,54,39,21,7,4,23,55,94,96,77,59,44,47,30,16,7,37,74]

mean(x(:))

x=x-mean(x)

Rx=xcorr(x)

y=x(end:-1:1)

subplot(2,1,1)stem(Rx)grid onxlabel('t')ylabel('Rx')title('xcorr()函数计算的自相关序列')

L=pow2(nextpow2(length(n1)+length(n2)-1))%确定FFT快速卷积的点数

Xk=fft(x,L)%计算x的L点FFT,结果为x(k)

Yk=fft(y,L)%计算y的L点FFT,结果为y(k)

Rk=Xk.*Yk%计算YK

Rn=ifft(Rk,L)%对RK调用IFFT,求得y1(n)

n=(n1(1)+n2(1)):(L+n1(1)+n2(1)-1)%确定n的自变量取值范围

subplot(2,1,2)stem(1:200,Rn(1:200))grid onxlabel('t')ylabel('R(n)')title('fft计算的自相关序列')

登录后复制

四、思考题

(1)对实验内容2中快速卷积基本算法的运算量进行分析,即乘法和加法运算次数。

(2)利用实例说明快速卷积基本算法的适用条件,即在什么情况下效率最高。

(3)实验内容3中,如何通过DFT的性质,减少一次ftt()函数的调用,提高计算效率?

五、总结

通过此次实验的练习,加深理解FFT 在实现数字滤波(或快速卷积)中的重要作用,更好的利用FFT进行数字信号处理,并且掌握了循环卷积和线性卷积两者之间的关系

时域卷颂让斗谈积与频域相乘是等效的;野销局

要用fft()算x和y的卷积的话就先补零算fft,再相乘,最后做ifft()就行了:

X = fft([x zeros(1,length(y)-1)])% x和y是要计算卷积的序列

Y = fft([y zeros(1,length(x)-1)])

z = ifft(X.*Y) % 结果存为z


欢迎分享,转载请注明来源:内存溢出

原文地址: http://outofmemory.cn/yw/12355365.html

(0)
打赏 微信扫一扫 微信扫一扫 支付宝扫一扫 支付宝扫一扫
上一篇 2023-05-24
下一篇 2023-05-24

发表评论

登录后才能评论

评论列表(0条)

保存