仲易+数字信号处理实验报告

1、北京理工大学20112012学年第一学期数字信号处理实验报告数字信号处理实验报告姓名：仲易班级：05110901学号：20091087实验1 基2-FFT算法实现一、实验目的1.掌握基2-FFT的原理及具体实现方法。2.编程实现基2-FFT算法。3.深刻理解FFT算法的特点。二、实验设备与环境计算机，matlab软件环境。三、 实验基础理论FFT是DFT的一种快速算法，能使DFT的计算大大简化，运算时间缩短。FFT利用了WNnk的三个固有特性。即对称性，周期性和可约性，将长序列的DFT分解为短序列的DFT，合并了DFT运算中的某些项，从而减少了DFT的运算量。FFT算法基本上可分为两大类，即按

2、时间抽取法和按频率抽取法。在实现FFT算法时，要重点考虑两个问题，注意数据的读取和存储：(1)输入输出的排序；（2）蝶形运算的实现。按时间抽取算法中输入反序输出顺序，按频率抽取算法中输入顺序输出反序；运算过程中的每一级都有N/2个蝶形运算构成，每一个蝶形运算单元中，两个节点变量运算后得到的结果为下一列相同位置的节点变量，而和其他节点变量无关，可以采用原位运算，节省存储单元。另外，蝶形运算中的复系数WNnk可以存储为能及时查阅的系数表，这样可以借阅运算量，但是需要N/2哥复数存储器。MATLAB中提供了用于计算FFT的函数fft，可将实验中所得到的结果与利用MATLAB中fft函数计算的结果相比

3、较，以此验证结果的正确性。4. 实验内容及实验过程1. 编程实现序列长度为N=8的按时间抽取的基2-FFT算法。给定一个8点序列，采用编写的程序计算其DFT，并与MATLAB中fft函数计算的结果相比较，以验证结果的正确性。代码如下：x=1 2 3 4 5 6 7 8;m=3; %X序列长度为2的3次幂N=8;nxd=bin2dec(fliplr(dec2bin(1:N-1,m)+1; %求倒序列y=x(nxd); for mm=1:m %对x做m次基2分解Nz=2mm;u=1; WN=exp(-i*2*pi/Nz); for j=1:Nz/2for k=j:Nz:Nkp=k+Nz/2; t=

4、y(kp)*u; %蝶形运算乘项y(kp)=y(k)-t; %加项y(k)=y(k)+t; endu=u*WN;endendyy1=fft(x) %对比matlab的ffty = Columns 1 through 3 36.0000 -4.0000 + 9.6569i -4.0000 + 4.0000i Columns 4 through 6 -4.0000 + 1.6569i -4.0000 -4.0000 - 1.6569i Columns 7 through 8 -4.0000 - 4.0000i -4.0000 - 9.6569iy1 = Columns 1 through 3 36

5、.0000 -4.0000 + 9.6569i -4.0000 + 4.0000i Columns 4 through 6 -4.0000 + 1.6569i -4.0000 -4.0000 - 1.6569i Columns 7 through 8 -4.0000 - 4.0000i -4.0000 - 9.6569i2、 编程实现序列长度为N=8的按频率抽取的基2-FFT算法。给定一个8点序列，采用编写的程序计算其DFT，并与MATLAB中fft函数计算的结果相比较，以验证结果的正确性。代码如下：L=2048;M=256;>> x=ones(1,L);>> n=0:

6、M-1;>> h=cos(0.2*pi*n);>> ticy1=conv(x,h);tocticX=fft(x,L+M -1);H=fft(h, L+M -1);Y=X.*H;y2=ifft(Y);toc运行结果如下：（分别为线形卷积和快速卷积）Elapsed time is 0.003352 seconds.Elapsed time is 0.007167 seconds.L=4096;M=256;>> x=ones(1,L);>> n=0:M-1;>> h=cos(0.2*pi*n);>> ticy1=conv(x,h

7、);tocticX=fft(x,L+M -1);H=fft(h, L+M -1);Y=X.*H;y2=ifft(Y);toc运行结果如下：（分别为线形卷积和快速卷积）Elapsed time is 0.003666 seconds.Elapsed time is 0.030134 seconds.3.将上述FFT程序推广大序列长度为N=2v的情况，要求利用原位运算。代码如下：L=2048;M=256;>> n=0:M-1;>> h=cos(0.2*pi*n);l=128;x=ones(1,l);q=0;ticH=fft(h,l+M-1);X=fft(x,l+M-1);Y

8、=X.*H;y=ifft(Y);yp=y;for i=2:16q=q+l;yq=zeros(1,q),y;yp=yp,zeros(1,l);yp=yq+yp;endtoc运行结果如下：Elapsed time is 0.018875 seconds.4.L=2048;M=256;n=0:M-1;h=cos(0.2*pi*n);x=ones(1,L);Lx=2048;M=256;M1=M-1;N=512;L=N-M1;tich=h,zeros(1,N-M); x=zeros(1,M1),x,zeros(1,N-1); a=floor (Lx+M1-1)/(L)+1; Y=zeros(1,N);

9、for k=0:a-1 xk=x(k*L+1:k*L+N); b=fft(xk,N); C=fft(h,N); Z=b.*C; Y(k+1,:)=ifft(Z,N); endY=Y(:,M:N)'Y=(Y(:)'toc运行结果如下：Elapsed time is 0.035329 seconds.实验2 利用DFT分析信号频谱一、 实验目的1、 加深对DFT原理的理解。2、 应用DFT分析信号的频谱。3、 深刻理解利用DFT分析信号频谱的原理，分析实现过程中出现的现象及解决方法。二、 实验设备与环境计算机、MATLAB软件环境三、 实验基础理论1、 DFT与DTFT的关系序列的

10、N点DFT实际上就是序列的DTFT在N个等间隔频率点上的采样样本。2、 利用DFT求DTFT方法一：利用差值运算 其中为内插函数方法二：在实际MATLAB计算中，采用增加数据的长度N，使得到的DFT谱线更加精细，用其包络近似计算DTFT。如果数据量不足，可采用补零来增加数据长度。3、 利用DFT分析连续时间信号的频谱1） 确定时域采样间隔T，得到离散序列2） 确定截取长度M，得到M点离散序列，这里为窗函数3） 确定频域采样点数N，要求4） 利用FFT计算离散序列的N点DFT，得到5） 由计算采样点的近似值。4、 可能用到的MATLAB函数与代码实验中DFT运算采用MATLAB中提供的函数fft

11、来实现。DTFT可以利用MATLAB矩阵运算的方法进行计算。四、 实验内容1、（1）n=0:3;x=2 -1 1 1;w=-pi:0.01*pi:pi;X=x*exp(-j*n'*w);subplot(211)plot(w,abs(X);axis tightsubplot(212)plot(w,angle(X)/pi);axis tight(2)在原基础上增加代码：n=0:3;x= 2 -1 1 1;N=4;k=0:3;Y=fft(x);w=0:0.01*pi:2*pi;X=x*exp(-j*n'*w);m=0:2*pi/4:3/2*pisubplot(211);plot(w,

12、abs(X);xlabel('Omega/pi');title('DTFT-Magnitude');axis tighthold onstem(m,abs(Y);subplot(212);plot(w,angle(X)/pi);xlabel('Omega/pi');title('DTFT-Phase');hold onstem(m,angle(Y)（3）在原本基础上设置代码subplot(211);Z=fft(x,64); % n=64DFTN=0:63;stem(abs(Z);subplot(212);stem(angle(Z)

13、2、考察序列（1），用DFT估计的频谱；将补零加长到长度为100点序列用DFT估计的频谱。要求画出相应波形。程序代码：程序输出结果：补零后：程序代码：程序输出结果：（2）时，用DFT估计的频谱，并画出波形。程序代码：程序输出结果：（3）根据实验结果，分析怎样提高频谱分辨率。 答：减小截取长度或减小采样间隔的方法可以提高频谱分辨率。3、已知信号其中=1Hz，=2Hz，=3Hz。从的表达式可以看出，它包含三个频率的正弦波，但是，从其时域波形来看，似乎是一个正弦信号，利用DFT做频谱分析，确定适合的参数，使得到的频谱的频率分辨率符合需要。程序代码： 程序输出结果：利用DFT近似分析连续时间信号的频谱

14、。分析采用不同的采样间隔和截取长度进行计算的结果，并最终确定适合的参数。程序代码：程序输出结果：通过改变步长可以达到改变采样间隔的结果，输出不同的频谱。实验5、脉冲响应不变法设计IIR数字滤波器一、 实验目的1、 掌握利用脉冲响应不变法设计IIR数字滤波器的原理及具体方法2、 加深理解数字滤波器和模拟滤波器之间的技术指标转化。3、 掌握脉冲相应不变法设计IIR数字滤波器的优缺点及适用范围。二、 实验设备与环境计算机、MATLAB软件环境三、 实验基础理论1、 基本原理：从时域响应出发，使数字滤波器的单位脉冲响应h（n）模仿模拟滤波器的单位冲激响应。2、 设计步骤：a. 确定数字滤波器性能指标。

15、b. 将数字滤波器频率指标转换成相应的模拟滤波器的频率指标c. 设计模拟滤波器实现方法：脉冲响应不变法 r,p,k=residue(b,a) b,a=residue(r,p,k)d. 模拟滤波器数字化实现方法：bz,az=impinvar(b,a,fs)四、 实验内容1、设采样频率fs=4khz,采用脉冲响应不变法设计一个三阶巴特沃斯数字低通滤波器，其3db截止频率为fc=1khz.首先设计巴特沃斯滤波器：N=3;Wc=0.25*pi;b,a=butter(N,Wc,'s')b = 0 0 0 0.4845a = 1.0000 1.5708 1.2337 0.4845因此，然后

16、脉冲响应不变法计算：bz,az=impinvar(b,a)bz = 1.1702 -1.4212 0.6078 -0.0945az = 1.0000 -1.6286 0.9286 -0.1879得到：画出模拟滤波器和数字滤波器的图像进行对比：H,w=freqs(b,a);subplot(221);plot(w/pi,abs(H);grid on;subplot(222);plot(w/pi,20*log10(abs(H)/max(abs(H);w=0:500*pi/500;H,w=freqz(bz,az);grid on;subplot(223);plot(w/pi,abs(H);grid o

17、n;subplot(224);plot(w/pi,20*log10(abs(H)/max(abs(H);grid on;2.设采样频率为fs=10khz，设计数字滤波器满足如下指标 通带截止频率：1khz，通带波动：1db， 阻带截止频率：1.5khz，阻带衰减：15db 分别采用巴特沃斯、切比雪夫I型、切比雪夫ii型和椭圆模拟原型滤波器及脉冲响应不变法进行设计。结合实验结果，分别讨论采用上述方法设计的数字滤波器是否都能满足给定指标要求，分析脉冲响应不变法设计IIR数字滤波器的优缺点及适用范围。（1）巴特沃斯滤波器为：Wp=0.1*pi;Ws=0.15*pi;Rp=1;As=15;N=ceil

18、(log10(10(Rp/10)-1)/(10(As/10)-1)/(2*log10(Wp/Ws)OC1=Wp/(10(Rp/10)-1)(1/(2*N)b,a=butter(N,OC1,'s')N = 6OC1 = 0.3516b = 0 0 0 0 0 0 0.0019a = 1.0000 1.3585 0.9227 0.3974 0.1141 0.0208 0.0019得到：bz,az=impinvar(b,a)bz = 1.0e-003 * 0.0000 0.0125 0.2573 0.5194 0.1637 0.0051 0az = 1.0000 -4.6517 9.

19、1374 -9.6829 5.8300 -1.8889 0.2570得到bz，az后，画出模拟滤波器和数字滤波器的图像进行对比：H,w=freqs(b,a);subplot(221);plot(w/pi,abs(H);grid on;subplot(222);plot(w/pi,20*log10(abs(H)/max(abs(H);w=0:500*pi/500;H,w=freqz(bz,az);grid on;subplot(223);plot(w/pi,abs(H);grid on;subplot(224);plot(w/pi,20*log10(abs(H)/max(abs(H);grid

20、on;H,w=freqs(b,a);subplot(221);plot(w/pi,abs(H);grid on;subplot(222);plot(w/pi,20*log10(abs(H)/max(abs(H);w=0:500*pi/500;H,w=freqz(bz,az);grid on;subplot(223);plot(w/pi,abs(H);grid on;subplot(224);plot(w/pi,20*log10(abs(H)/max(abs(H);grid on;（2）切比雪夫型>> Wp=0.2*pi;>> Ws=0.3*pi;>> Rp=

21、1;>> As=15;>> e=>> e=(10(Rp/10)-1)(1/2)e =0.5088>> A=10(As/20)A = 5.6234>> N=ceil(acosh(A*A-1)(1/2)/e)/acosh(Ws/Wp)N = 4>> b,a=cheby1(N,Rp,Wp,'s')b = 0 0 0 0 0.0383a =1.0000 0.5987 0.5740 0.1842 0.0430（3）切比雪夫型：>> Wp=0.2*pi;>> Ws=0.3*pi;>>

22、 Rp=1;>> As=15;>> N=4;>> b,a=cheby2(N,As,Ws,'s')b =0.1778 -0.0000 1.2637 -0.0000 1.1225a =1.0000 2.3696 3.0322 1.9925 1.1225（4）椭圆模拟原型滤波器：>> Wp=0.2*pi;>> Ws=0.3*pi;>> Rp=1;>> As=15;>> e=0.5088;>> A=5.6234;>> k1=e/(A*A-1)(1/2)k1 =0.09

23、19>> k=Wp/Wsk = 0.6667>> K,E=ellipke(k)K =2.0290E =1.2612>> K1,E1=ellipke(1-k1*k1)(1/2)K1 =4.1217E1 =1.0077>> K2,E2=ellipke(k1)K2 = 1.6089E2 =1.5340>> K3,E3=ellipke(1-k*k)(1/2)K3 =2.1483E3 = 1.2140>> N=(K*K1)/(K2*K3)N = 2.4195>> N2=2;>> b,a=ellip(N2,Rp

24、,As,Ws,'s')b = 0.1778 0 0.9379a =1.0000 0.9120 1.0523结果分析：（1） 当给定的设计指标相同时，选用椭圆滤波器所要求的结束N最低，切比雪夫滤波器 次之，巴特沃斯滤波器最高；（2） 从通带的相位相应来看，椭圆滤波器虽然提供了最优秀的幅度平方相应，但通带上 的相位相应非线性较大，而巴特沃斯滤波器在通带上具有相当的线性相位，切比雪夫滤波器的相位特征介于两者之间。实验七、窗函数法设计FIR数字滤波器一、实验目的掌握窗函数法设计FIR数字滤波器的原理及具体方法。二、实验设备及环境计算机、matlab软件及环境三、实验基础理论1.基本原理

25、 窗函数设计法的基本思想为，首先选择一个适当的理想的滤波器Hd（ejw），然后用窗函数截取它的单位脉冲响应hd（n），得到线性相位和因果的FIR滤波器。这种方法的重点是选择一个合适的窗函数和理想滤波器，使设计的滤波器的单位响应尽力逼向滤波器的单位响应。2.设计步骤（1）给定理想滤波器的频率响应Hd（ejw），再通带上具有单位增益和线性相位，在阻带上具有零响应。（2）确定这个滤波器的单位脉冲响应hd=sincn-n-为了得到一个h(n)长度为N的因果相位FIR滤波器，我们令=N-12（3）用窗函数截取hd（n）得到所设计FIR数字滤波器h(n)h(n)= hd（n）w（n）四、实验内容1.设计一

26、个数字低通滤波器，其技术指标如下 p=0.2, Rp=0.25db st=0.3, As=50db分别采用矩形窗、汉宁窗、海明窗、布莱克曼窗、凯瑟窗设计该滤波器。结合实验结果，分别讨论采用上述方法设计的数字滤波器是否都能满足给定指标要求。2. 设计一个数字低通滤波器，其技术指标如下 下阻带边缘：st1=0.2, As=60db 下通带边缘：p1=0.35, Rp=1db 上通带边缘：p2=0.65, Rp=1db 上阻带边缘：st2=0.8, As=60db五、代码确定长度：wp=0.2*pi;>> wst=0.3*pi;>> width=wst-wp;>>

27、 N=ceil(1.8*pi/width)+1N = 19确定hd(n)>> n=0:(N-1);>> wc=(wp+wst)/2;>> alpha=(N-1)/2;>> hd=(wc/pi)*sinc(wc/pi)*(n-alpha)>> wboxcar=boxcar(N)'>> h=hd.*wboxcarsubplot(221);>> stem(n,hd,'filled');axis tight;xlabel('n');ylabel('hd(n)');

28、>> Hr,w1=zerophase(h);>> subplot(222);plot(w1/pi,Hr);>> axis tight;xlabel('omegapi');ylabel('omega');>> subplot(223);stem(n,h,'filled');axis tight;xlabel('n');ylabel('h(n)');>> H,w=freqz(h,1);>> subplot(224);>> plot(w/

29、pi,20*log10(abs(H)/max(abs(H);>> xlabel('omegapi');ylabel('dB')>> grid on由上图可清楚看出：加矩形窗不满足指标用海明窗：wp=0.2*pi;wst=0.3*pi;width=wst-wp;N=ceil(6.6*pi/width)+1N = 67n=0:(N-1);wc=(wp+wst)/2;alpha=(N-1)/2;hd=(wc/pi)*sinc(wc/pi)*(n-alpha);whamming=hamming(N)'h=hd.*whamming;>

30、> subplot(221);stem(n,hd,'filled');axis tight;xlabel('n');ylabel('hd(n)');Hr,w1=zerophase(h);subplot(222);plot(w1/pi,Hr);axis tight;xlabel('omegapi');ylabel('omega'); subplot(223);stem(n,h,'filled');axis tight;xlabel('n');ylabel('h(n)

31、9;);H,w=freqz(h,1);subplot(224);plot(w/pi,20*log10(abs(H)/max(abs(H);xlabel('omegapi');ylabel('dB')grid on由上图可清楚看出：加海明窗满足指标 汉宁窗：代码：wp = 0.2*pi;wst = 0.3*pi; tr_width = wst-wp; N = ceil(6.2*pi/tr_width)+1; n = 0:(N-1);wc = (wp+wst)/2; alpha = (N-1)/2; hd = (wc/pi)*sinc(wc/pi)*(n-alpha

32、); w_hanning = hanning(N)' h = hd.*w_hanning; subplot(221);stem(n,hd,'filled'); axis tight;xlabel('n');ylabel('hd(n)'); Hr,w1 = zerophase(h); subplot(222);plot(w1/pi,Hr); axis tight;xlabel('omega/pi');ylabel('H(omega)'); subplot(223);stem(n,h,'filled'); axis tight;xlabel('n');ylabel('h(n)'