用DFT计算线性卷积

实验一. 用DFT 计算线性卷积
1． 实验要求
1.1令w(n)为一正弦加白噪声信号，长度为500，h(n)是一个低通FIR 滤波器，试用DFT 与叠接相加法计算线性卷积。

h(n)使用fir1.m 设计：
h=fir1(10,0.3,hanning(11)); %设计低通滤波器，得到h(n);
1.2设计一个按照时间抽取的基2快速傅里叶变换（基2FFT-DIT ）。

2． 实验原理
w(n)为长序列，将其分解为若干个子序列，利用叠接相加法进行卷积。

设x(n)为一个M 点序列，h(n)为一个L 点序列，y(n)=x(n)*h(n),即y(n)是x(n)和h(n)的线性卷积，那么y(n)是一(M+L-1)点的序列。

又由于DFT 对应循环卷积而不是线性卷积，能否用DFT 来计算两个序列的线性卷积呢？答案是肯定的。

由DFT 的性质有：
()()[()()]x n h n IDFT X k H k ⊗=
式中x(n),y(n)都是N 点的序列，循环卷积的结果也是N 点序列，显然X(k),H(k)也是N 点序列，现希望
()()()[()()][()]y n x n h n IDFT X k H k IDFT Y k =*==
因为y(h)是（M+L-1）点序列，因此，Y(k)也必须是（M+L-1）点序列，相应的X(k),H(k)也都应当是M+L-1点序列，而且X(K),H(K)对应的时域序列x(n)，h(n)也必须是M+L-1点的序列。

只有这样，由Y(k)作逆变换所得到的y(n)才能保证x(n)和h(n)的线性卷积，具体步骤如下：
1) 对M 点序列x(n)及L 点序列h(n)分别作扩展，构成新序列()x n '，()h n '，他们的长度都是M+L-1点，即
()0,1,,1()0
,2x n n M x n n M M L =-⎧'=⎨=+-⎩ ()0,1,,1()0
,2h n n L h n n L M L =-⎧'=⎨=+-⎩
2) 认为()x n '，()h n '各是周期序列()x n '，()h n '的一个周期，周期长度为M+L-1,可得
()()(),
1y n x n h n N M L '''=⊗=+-
而
()()()*()y n y n x n h n '==
3) 若用DFT 求y(n),则
()()[()()]y n y n IDFT X k H k '''==
式中()X n '，()H n '分别是()x n '，()h n '的DFT 。

对于长序列的卷积，可以把长序列分成短序列分别作卷积，再按一定的规则首位相加。

相应的方法有叠接相加法和叠接舍去法。

3． 实现代码及结果
h=fir1(10,0.3,hanning(11));
N=500;p=0.05;f=1/16;
u=randn(1,N)*sqrt(p);
s=sin(2*pi*f*[0:N-1]);
x=u(1:N)+s;
y=fftfilt(h,x);
subplot(211)
plot(x);
subplot(212)
plot(y);
clc;close all;clear;format compact;
%输入数据并计算常量
xn=[0,1,2,3,4,5,6,7];%可取任意序列
M=nextpow2(length(xn)), N=2^M,
for m=0:N/2-1;%旋转因子指数范围
WN(m+1)=exp(-j*2*pi/N)^m;%计算旋转因子end
A=[xn,zeros(1,N-length(xn))]; %数据输入disp('输入到各存储单元的数据:'),disp(A); %数据倒序操作
J=0;%给倒序数赋初值
for I=0:N-1;%按序交换数据和算倒序数
if I<J;%条件判断及数据交换
T=A(I+1);A(I+1)=A(J+1);A(J+1)=T;
end
%算下一个倒序数
K=N/2;
while J>=K;
J=J-K;K=K/2;
end
J=J+K;
end
disp('倒序后各存储单元的数据:'),disp(A);
%分级按序依次进行蝶形运算
for L=1:M;%分级计算
disp('运算级次:'),disp(L);
B=2^(L-1);
for R=0:B-1;%各级按序蝶算
P=2^(M-L)*R;
for K=R:2^L:N-2;%每序依次计算
T=A(K+1)+A(K+B+1)*WN(P+1);
A(K+B+1)=A(K+1)-A(K+B+1)*WN(P+1); A(K+1)=T;
end
end
disp('本级运算后各存储单元的数据:'),disp(A); end
disp('输出各存储单元的数据:'),Xk=A,
disp('调用fft函数运算的结果:'),fftxn=fft(xn,N),。