用DFT计算线性卷积

用DFT计算线性卷积

1 基本原理

1.1用 DFT实现线性卷积的原理

根据线性卷积的原理 :

x[n]  y[n] D F T X(e jω)Y( e jω),

且， w[n] = x[n]  y[n] 可用下式求得：

－1

F

{X(e j)Y( e j)}

ωω线性与圆周卷积分别由下式给出

w[n]x[n]y[n]mx[m]y[nm]m0wp[n]x[n] y[n]N x[m]y[((nm))N1N]RN[n]其中 x[n] : 0≤n ≤P －1  0≤m ≤P －1 y[n]: 0≤n ≤L －1  0≤n － m ≤L －1 w[n]的最大长度为 ：L+P-1，单 wp[n] 的长度为 N。 当N≥L+P －1 , wp[n] = w[n]; 当 N ≤ L+P －1, wp[n] ≠ w[n];

所以要使圆周卷积等于线性卷积而不产生混叠的必要条件为：N ≥ P+L－1 即线性与圆周卷积一致的样本为: P+L － N－1≤ n ≤N －1

1.2 重叠保留法原理

设h(n)的点数为M，信号x(n)为很长的序列。我们将x(n)分解为很多段，每段为L点，L选择成和M的数量级相同，用xi(n)表示x(n)的第i段：

要求x(n)和h(n)的卷积时，若x(n)的点数很多，远大于h(n)的点数M时，通常不允许等x(n)全部采集齐后再进行卷积，否则，使输出相对于输入有较长的延时。因此需要采用

2

分段卷积或称分段过滤的办法，即将x(n)分成点数和h(n)相仿的段，分别求出每段的卷积结果，然后用一定方式把它们合在一起，便得到总的输出，一种分段卷积的方法就是重叠保留法。

设h(n)的点数为M，信号x(n)为很长的序列。我们将x(n)分解为很多段，每段为L点，L选择成和M的数量级相同，取N=L+M-1用xi(n)表示x(n)的第i段。重叠保留法先将x(n)分段，每段L=N-M+1个点，由于xi(n)*h(n)为L+M-1 点，故先对xi(n)及h(n)补零值点，补到N点不同之处是，序列中补零处不补零，而在每一段的前边补上前一段保留下来的（M-1）个输入序列值， 组成L+M-1点序列xi(n)，如图8-27（a）所示。如果L+M-1<2m， 则可在每段序列末端补零值点，补到长度为2m，这时如果用

DFT实现h(n)和xi(n)圆周卷积，

则其每段圆周卷积结果的前（M-1）个点的值不等于线性卷积值，必须舍去。

为了说明以上说法的正确性，我们来看一看图8-28。任一段xi(n)（为N点）与h(n)（原为M点，补零值后也为N点）的N点圆周卷积

yi(n)xi(n)xi(m)h((nm))NRN(n)N h(n)m0N1由于h(m)为M点，补零后作N点圆周移位时，在n=0,1,…,M-2的每一种情况下，h((n-m))NRN(m)在0≤m≤N-1范围的末端出现非零值， 而此处xi(m)是有数值存在的，图8-28（c），（d）为n=0, n=M-2的情况，所以在0≤n≤M-2

这一部分的yi(n)值中将混入xi(m)尾部与h((n-m))N·RN(m)尾部的乘积值，从而使这些点的yi(n)不同于线性卷积结果。但是从n=M-1开始到n=N-1，h((n-m))NRN(m)=h(n-m（)如图8-28（e），（f）所示），圆周卷积值完全与线性卷积值一样，yi(n)就是正确的线性卷积值。因而必须把每一段圆周卷积结果的前(M-1)个值去掉， 如图8-28（g）所示。

因此，为了不造成输出信号的遗漏，对输入分段时，就需要使相邻两段有M-1个点重叠（对于第一段，即x0(n)，由于没有前一段保留信号，则需要在序列前填充M-1 个零值点），这样，若原输入序列为x′(n)（n≥0 时有值），则应重新定义输入序列

x[ni(NM1)]xi(n)00≤n≤M-2 M-1≤n

00≤n≤N-1

x(n)x ' [ n   ( M  1 )] 其他n （ i=0, 1, … ）

3

在这一公式中，已经把每一段的时间原点放在该段的起始点，而不是x(n)的原点。这种分段方法示于图8-27中，每段xi(n)和h(n)的圆周卷积结果以yi(n)表示，如图 8-27（b）所示，图中已标出每一段输出段开始的(M-1)个点，0≤n≤M-2部分舍掉不用。把相邻各输出段留下的序列衔接起来，就构成了最后的正确输出， 即 式中:

y(n)yi'[ni(NM1)]i0 '  y i ( n ) M-1≤n≤N-1

yi(n)  0 其他n

(c)0这时，每段输出的时间原点放在yi ′ (n)的起始点，而不是y(n)的原点。

h(m)h((M－2－m))NRN(m)(a)0M－1xi(m)N－1m(d)0n＝M－2N－1h((M－1－m))NRN(m)m(e)0n＝M－1N－1h((N－1－m))NRN(m)m(b)0h((0－m))NRN(m)N－1m(f )n＝0N－1m0yi(n)n＝N－1N－1m(g)0M－2N－1n图1-1 用保留信号代替补零后的局部混叠现象

x0(n)M－2N－1ny0(n)M－2LM－10N－1ny1(n)N－1N－1n

(a)0x1(n)0M－2L0(b)M－2M－1y2(n)M－2nx2(n)0M－2LN－1nM－1N－1n图1-2 重叠保留法示意图

4

2函数说明

由于是用重叠保留法计算线性卷积，根据重叠保留法的原理，在MATLAB软件的环境下可以创造一个相应的函数求x(n)和h(n)的线性卷积。设h(n)的点数为M，信号x(n)为很长的序列。把x(n)分成长为L 的小段xi(n),每段与前一段重叠M-1 个样本的多段,保留后面的(L-M+1)个输出样本,最后把这些输出连成一个序列，结果即为x(n)和h(n)的线性卷积。

下面的函数程序给出这种重叠保留分段卷积法的算法 (MATLAB 代码)： %dft.m

function [Xk]=dft(xn,N) %计算离散傅里叶变换 %----------------------------------- %[Xk]=dft(xn,N)

%Xk=在0<=n<=N-1间的DFT系数数组 %xn=N点有限持续时间序列 %N=DFT的长度 %

n=[0:1:N-1]; %n的行向量 k=[0:1:N-1]; %k的行向量 WN=exp(-j*2*pi/N); %Wn的因子

nk=n'*k; %产生一个含nk值的N乘N维矩阵 WNnk=WN.^nk; %DFT矩阵 Xk=xn*WNnk; %DFT系数的行向量 %idft.m

function [xn]=idft(Xk,N) %计算逆离散傅里叶变换 %----------------------------- %[xn]=idft(Xk,N)

5

%xn=在0<=n<=N-1间的N点有限持续时间序列 %Xk=在0<=k<=N-1间的DFT系数数组 %N=DFT的长度 %

n=[0:1:N-1]; %n的行向量 k=[0:1:N-1]; %k的行向量 WN=exp(-j*2*pi/N); %Wn的因子

nk=n'*k; %产生一个含nk值的N乘N维矩阵 WNnk=WN.^(-nk); %IDFT矩阵 xn=(Xk*WNnk)/N; %IDFT的行向量

%circonvf.m

function y=circonvf(x1,x2,N)

%在x1 和x2：（频域）之间的N 点圆周卷积 %----------------------------------------- %[y]=circonvf(x1,x2,N) %y=包含圆周卷积的输出序列 %x1=长度N1<=N的输入序列 %x2=长度N2<=N的输入序列 %N=循环缓冲器的大小 %方法y(n)=IDFT(X1(k).X2(k)) %检查x1的长度： if length(x1)>N

error('N必须>=x1的长度') end

%检查x2的长度： if length(x2)>N

error('N必须>=x2的长度')

6

end

x1=[x1 zeros(1,N-length(x1))]; x2=[x2 zeros(1,N-length(x2))]; X1=dft(x1,N); X2=dft(x2,N); X=X1.*X2; Y=idft(X,N) y=real(Y);

%ovrlpsav.m

function[y]=ovrlpsav(x,h,N) %用重叠保留法作分段卷积 %[y]=ovrlpsav(x,h,N) %y=输出序列 %x=输入序列 %h=脉冲响应 %N=段长

Lenx=length(x);M=length(h); M1=M-1;L=N-M1; h=[h zeros(1,N-M)]; x=[zeros(1,M1),x,zeros(1,N-1)]; K=floor((Lenx+M1-1)/(L)); Y=zeros(K+1,N);

for k=0:K xk=x(k*L+1:k*L+N); Y(k+1,:)=circonvt(xk,h,N); end

Y=Y(:,M:N)'; y=(Y(:))';

%输入数据x的长度及脉冲响应h的长度 %各段搭接长度M1，有效数据长度L %将h延长至循环长度N %把x前面加上(M-1)个零 %段数 %各段进行卷积 %Y中各行去掉前(M-1)个样本 %转成输出

7

3 设计程序

1,0,1,0,1用重叠保留法计算线性要求对序列x(n)=n+1（0n34）和h(n)卷积，利用MATLAB软件编程绘制出每一段卷积结果图和 x(n)*h(n)图。在已有函数[y]=ovrlpsav(x,h,N)的情况下，根据要求编写程序如下：

4 仿真结果 5 心得体会 6 参考文献

[1] 余成波等编，数字信号处理及MATLAB实现（第二版），北京：清华大学出版社，2008 [2] 陈怀琛编著，数字信号处理教程——MATLAB释义与实现，北京：电子工业出版社2004 [3] 万永革编著，数字信号处理的MATLAB实现，北京：电子工业出版社2007 [4] 薛年喜编著，MATLAB在数字信号处理中的应用，北京：清华大学出版社，2008 [5] 李正周编著，MATLAB数字信号处理与应用，北京：清华大学出版社，2008

8