实验四 利用FFT 实现快速卷积
一、实验要求
加深理解FFT 在实现数字滤波(或快速卷积)中的重要作用,更好的利用FFT 进行数 字信号处理。掌握循环卷积和线性卷积两者之间的关系。
二、实验原理
数字滤波器根据系统的单位脉冲响应h(n) 是有限长还是无限长可分为有限长单位脉冲响应(Finite Impulse Response)系统(简记为FIR系统)和无限长单位脉冲响应(Infinite Impulse Response)系统(简记为IIR系统)。
对于FIR滤波器来说,除了可以通过数字网络来实现外,也可以通过FFT的变换来实现。首先我们知道,一个信号序列x(n) 通过FIR滤波器时,其输出应该是x(n)与h(n)的卷积:
当h(n)是一个有限长序列,即h(n)是FIR滤波器,且时
在数字网络类的FIR滤波器中,普遍使用的横截型结构就是按这个卷积公式构成的。
应用FFT实现数字滤波器实际上就是用FFT来快速计算有限长度序列的线性卷积。
这种方法就是先将输入信号x(n) 通过FFT变换为它的频谱采样值X(k),然后和FIR滤波器的频响采样值H(k)相乘,H(k)可事先存放在存储器中,最后将乘积H(k)X(k)通过快速傅里叶反变换(简称IFFT)还原为时域序列,即得到输出y(n)。
现以FFT求有限长序列间的卷积,和求有限长度序列与较长序列间的卷积为例来讨论FFT的快速卷积方法。
a.序列x(n)和h(n)的长差不多。设x(n)的长为N1,h(n)的长为N2,要求
用FFT完成这一卷积的具体步骤如下:
①为使两有限长序列的线性卷积可用其循环卷积代替而不发生混叠,必须选择循环卷积长度,若采用基2-FFT完成卷积运算,要求(m为整数)。
②用补零方法使x(n)和h(n)变成列长为N的序列。
③用FFT计算x(n)和h(n)的N点离散傅里叶变换
④完成X(k)和H(k)乘积,
⑤用FFT计算Y(k)的离散傅里叶反变换得
b.当x(n)长度很长时,即,通常不允许等x(n)全部采集齐后再进行卷积,否则使输出相对于输入有较长的延时,另外,若太大,h(n)要补上太多的零点,很不经济,且FFT的计算时间也要很长。为此,采用分段卷积的方法,即把x(n)分成长度与h(n)相仿的小段,分别求出每段卷积的结果,然后用相应的方式把它们结合起来,便是总的输出。分段卷积方法主要有两种,即重叠相加法和重叠保留法。具体内容请参考数字信号处理教材中“快速离散傅里叶变换”一章中的线性卷积的FFT算法部分,本实验这部分不作重点要求。
三、主要实验仪器及材料
微型计算机、Matlab6.5 教学版、TC 编程环境。
四、实验内容
(1)用Matlab 或C语言编制信号产生子程序,产生典型信号供谱分析用;
(2)对给出信号逐个进行谱分析,绘出序列和幅频特性曲线;
(3)设计利用快速傅里叶变换FFT 计算线性卷积的程序;
(4)对结果进行分析;
(5)完成实验报告。
五、实验结果
1.计算序列的线性卷积。
1)在命令行窗口输入如下程序
x=[2, 1, 3, 2, 1, 5, 1]; h=[1, 2, -1, -3];
N=length(x) + length(h) - 1;
n=0:N-1;
x=[x, zeros(1, N-length(x))];
h=[h, zeros(1, N-length(h))];
X=fft(x); H=fft(h);
Y=X.*H; y=ifft(Y);
subplot(221), stem(n, x, '.'); xlabel('n'); ylabel('x(n)'); title('序列x(n)'); grid;
subplot(222), stem(n, h, '.'); xlabel('n'); ylabel('h(n)'); title('序列h(n)'); grid;
subplot(223), stem(n, y, '.'); xlabel('n'); ylabel('y(n)'); title('序列y(n)'); grid;
2)运行后得如下图
2.已知,用重叠保留法(采用FFT计算)。求,该分段长度N=6。
1)编写hsolpsav.m文件
function[y]=hsolpsav(x, h, N)
N=2^(ceil(log10(N)/log10(2)));
Lx=length(x); M=length(h);
M1=M-1; L=N-M1;
H=fft(h, N);
K=ceil(Lx/L);
x=[zeros(1, M1), x, zeros(1, N-1)];
y=zeros(K, N);
for k=0:K-1
Xk=fft(x(k*L+1:k*L+N));
y(k+1, :)=real(ifft(Xk.*H));
end
y=y(:, M:N)'; y=(y(:))';
2)在命令行窗口输入
n=[0:12]; x=13-n; h=[2, -1, 1]; N=6;
y=hsolpsav(x, h, N);
3)运行后得到如下结果
3.已知,求4点圆周卷积。
1)编写circonvt.m文件
function y = circonvt(x1, x2, N)
if length(x1) > N
error('N必须 >= x1的长度')
end
if length(x2) > N
error('N必须 >= x2的长度')
end
x1=[x1 zeros(1, N-length(x1))];
x2=[x2 zeros(1, N-length(x2))];
X1=fft(x1, N);
X2=fft(x2, N);
X=X1.*X2;
y=ifft(X, N);
y=real(y);
2)在命令行窗口输入
x1=[2 4 3 1], x2=[2 1 3];
y=circonvt(x1, x2, 4)
3)得到如下结果
六、实验总结
通过本次实验我学会了FFT的分段卷积法——重叠相加法实现线性卷积的MATLAB程序,时域作重叠保留法实现线性卷积的MATLAB程序以及圆周卷积的频域方法;我的把数学算法转换为程序的能力得到提高;今后编程能力有待提高。