圆周卷积的积分算法
2 圆周卷积原理
2.1 圆周卷积定义
有两个有限长序列x 1(n ) 和x 2(n ) 。设x 1(n ) 和x 2(n ) 的长度分别为N1和N2,N=max[N1, N2]。
⎧x (n ) 0≤n ≤N 1-1
x 1(n ) =⎨1 N 1≤n ≤N -1⎩0 X 3(k )=X 1(k )·X 2(k ) 则
⎧x (n ) 0≤n ≤N 2-1
x 2(n ) =⎨2
N 2≤n ≤N -1⎩0
若x 1(n ) 、x 2(n ) 的离散付里叶变换分别为X 1(k )、X 2(k ),且有
x 3(n )=IDFT[X 3(k )] =[∑x 2(m )x 1((n -m ))N ]R N (n )
m =0
上式即为序列x 1(n ) 与x 2(n ) 的圆周卷积,习惯表示为 x 3(n )=x 1(n ) ⊙x 2(n ) 。 可以看出,圆周卷积和周期卷积过程一样,只不过这里要取结果的主值序列。公式中的x 1((n -m ))N 只在m=0到N-1范围内取值,因而它就是圆周移位,所以这一卷积和称为圆周卷积。圆周卷积与周期卷积之间的关系,就是有限长序列圆周卷积结果的周期延拓,等于它们周期延拓后的周期卷积。换句话说,周期卷积的主值序列,是各周期序列主值序列的圆周卷积。周期卷积得到是周期序列,圆周卷积得到的是有限长序列,而且长度等于参加卷积的序列的长度。
N -1
2.2 圆周卷积计算过程
具体步骤如下:
(1) 在二元坐标上做出x 1(m ) 与x 2(m ) ; (2) 把x 2(m ) 沿着纵坐标翻转,得到x 2(-m ) ; (3) 对x 2(-m ) 做圆周移位,得到x 2((n -m ))N R N (n );
(4) x 1(m ) 与x 2((n -m ))N R N (n )对应相同的m 的值进行相乘,并把结果相加,得到对应于
自变量n 的一个x 3(m ) ;
(5) 换另一个n, 重复以上两步,直到n 取遍0到N-1所有的值,得到完整的x 3(n ) 。
3 重叠相加法原理
实际应用中通常参加卷积的两个序列的长度相差较大,这样长度较小的就需要补很多0,就需要较大的存储空间,运算时间也会变长,常用的办法有重叠相加法和重叠保留法。本次设计采用重叠相加法。
设x(n)为无限长序列,并假定x i (n)表示第i 段x (n )序列:
⎧x (n ) iN 2≤n ≤(i +1) N 2-1
x i (n ) =⎨
0⎩
则输入序列可表为:
x (n ) = ∞i=−∞x i n
∞
于是输出可分解为:y (n )=x (n ) *h (n ) = ∞i=−∞x i n ∗h n = i=−∞y i n
其中y i (n )=x i (n ) *h (n ) 由此表明,只要将x (n )的每一段分别与h (n )卷积,然后再将这些卷积结果相加起来就可得到输出序列,这样,每一段的卷积都可用上面讨论的快速卷积来计算。先对h (n )及x i (n )补零,补到具有N 点长度,N=N1+N2-1。 一般选择N=2M ,然后用基2 FFT算法通过正反变换计算
y i (n )=xi (n )*h(n )
由于y (的长度为N ,而x (的长度为N 2,因此相邻两y (序列必然有N-N 2=N1-1i n )i n )i n )点发生重叠,这个重叠部分应该相加起来才能构成最后的输出序列。 计算步骤:
a. 事先准备好滤波器参数H (k )=DFT[h(n )],N 点 b. 用N 点FFT 计算X i (k )=DFT[xi (n)]
c.Y i (k)=Xi (k)H(k)
d. 用N 点IFFT 求 yi (n )=IDFT[Yi (k)] e. 将重叠部分相加
图3.1重叠相加法的分段示意图
4 MATLAB 程序设计
基于以上对原理的分析,这里采用FFT 算法实现重叠相加法的圆周卷积计算。程序主要由计算圆周卷积的函数以及重叠相加函数和主程序三部分组成。
4.1 程序设计思路
函数conv (x1,x2,L)设计
(1) x1(n)进行N 点快速傅里叶变换得X1k (2) x2(n)进行N 点快速傅里叶变换得X2k (3) 进行频域相乘Yk=X1k*X2k
(4) 对Yk 进行反变换得到时域卷积结果y(n) 函数overlap_add(x,h,N)设计
(1)首先取圆周卷积的周期L (即进行L 点的快速傅里叶变换) (2)计算每一分段的大小N
(3)填充序列使得循环中对序列的索引不会超出范围 (4)计算分段数K
(5)对序列进行分段调用conv ()函数计算圆周卷积 (6)各段重叠相加 (7)取出实际的输出序列
4.2 圆周卷积部分
计算圆周卷积函数的流程图:
程序:
function y = circular_conv( x1,x2,L ) %利用循环卷积计算线性卷积
%循环卷积采用频域计算方法,以FFT 代替DFT ,降低运算量 x1k= fft(x1,L);%x1做L 点FFT x2k= fft(x2,L);%x2做L 点FFT yk = x1k.*x2k;%频域相乘
y = ifft(yk);%FFT反变换得循环卷积结果 end
4.3 重叠相加部分
重叠相加法程序流程图:
程序:
function y = overlap_add(x1,x2,N) %重叠相加法实现
%将高点数DFT 转化为低点数DFT M = length(x2);%获得x2(n)的长度 if N
T = ceil(Lx/N);%确定分段数
t = zeros(1,M-1);%生成1行M-1列的全零矩阵
x1 = [x1,zeros(1,(T+1)*N-Lx)];%不足部分补0 y = zeros(1,(T+1)*N);%定义输出序列,初始化为全0 for i=0:1:T xi = i*N+1;
x_seg = x1(xi:xi+N-1);
y_seg = circular_conv(x_seg,x2,L);%计算每段的圆周卷积 y_seg(1:M-1) = y_seg(1:M-1)+t(1:M-1);%重叠相加 t(1:M-1) = y_seg(N+1:L); y(xi:xi+N-1) = y_seg(1:N); end
y=y(1:Lx+M-1); end
4.4 主程序
x1=[1:1:13]; x2=[-2,1,3]; L=15;
m=overlap_add(x1,x2,L);%重叠相加法得到的结果 n=circular_conv(x1,x2,L);%线性卷积直接得到的结果 subplot(4,1,1); a = 1:13; stem(a,x1) title('x1序列')
subplot(4,1,2); b = 1:3; stem(b,x2) title('x2序列')
subplot(4,1,3);
c = 1:15; stem(c,m)
title('基于重叠相加法的圆周卷积结果')
5 MATLAB 实现
图5.1 圆周卷积m 文件
图5.2 重叠相加法m 文件
图5.3 主程序m 文件
图5.4 程序运行结果
图5.5 基于重叠相加法的圆周卷积结果表格
心得体会
数字信号处理是我们电信专业非常重要的一门课程,学好它对以后的进一步学习或者工作都有至关重要的作用,经过了一个学期的理论课程的学习,我感觉自己还停留在做题目上,而对理论知识的应用却知之甚少,因此,这次课程设计对我来说是一个非常好的机会来进一步认识理论知识的具体意义,怀着这种心态,我从对理论知识的回顾到设计仿真,每一步都非常认真的对待,我从中也学到了不少。
这次课程设计主要用到了matlab 软件,这款软件在对数字信号处理的建模、编程、分析、实现等方面功能非常强大,由于之前在实验课中我们曾经接触过matlab ,因此这次做起来相对轻松一些。我选做的题目是基于重叠相加法的圆周卷积,在开始进行程序设计之前,我先对理论知识进行了回顾,在熟悉了圆周卷积以及重叠相加法的原理以及matlab 的应用环境后,我开始了软件的设计,虽说之前对matlab 有所了解,但对软件里面的程序包并不了解,所以编程的过程中也遇到过很多问题,比如找不到需要用的函数,但通过查看软件自带的帮助信息我很快能够找到自己想要的东西,在熟悉了原理后,程序设计并不复杂,但程序的调试却花了很长时间,其中也遇到了一些莫名其妙的问题,调试成功后第二次再打开就不行了,经过了长时间的调试,并且在老师和同学的帮助下,我换了个软件终于解决了。
虽然调试的过程并不顺利,但成功后的喜悦是无穷的。我从中学到了很多东西,这些都是书本上很难遇到的,当然最终能够获得成功也离不开老师的帮助,所以在此感谢老师对我的帮助
10
致谢
在这次课程设计过程中,我得到了许多人的帮助。
首先我要感谢吴巍老师在课程设计上给予我的指导、提供给我的支持和帮助,这是我能顺利完成这次报告的主要原因,更重要的是老师帮我解决了许多技术上的难题,让我能做得更加完善。在此期间,我不仅学到了许多新的知识,而且也开阔了视野,提高了自己的设计能力。
其次,我要感谢帮助过我的同学,同时也感谢学院为我提供良好的做课程设计的环境。 最后再一次感谢所有在设计中曾经帮助过我的良师益友和同学!
11
参考文献
[1]刘泉, 阙大顺, 郭志强. 数字信号处理. 北京:电子工业出版社,2009
[2]唐昌建.Matlab 编程基础与应用. 四川:四川大学网络教育学院,2003
[3]陈怀琛. 数字信号处理教程-Matlab 释疑与实现. 北京:电子工业出版社,2004
[4]张志涌.Matlab 教程. 北京:北京航空航天大学出版社,2005
[5]曹弋.matlab 教程与实例. 北京:机械工业出版社,2008
12