数字信号处理实践报告

《数字信号处理实践报告》题目数字信号处理实践报告专业电子信息工程学生姓名准考证号指导教师2012年5月一、目的与要求掌握《数字信号处理基础》课程的基本理论；掌握应用MATLAB进行数字信号处理的程序设计方法。二、主要内容设计题目及设计要求：已知低通数字滤波器的性能指标如下：，，，选择合适的窗函数，设计满足上述指标的数字线性相位FIR低通滤波器。用一个图形窗口，包括四个子图，分析显示滤波器的单位冲激响应、相频响应、幅频响应和以dB为纵坐标的幅频响应曲线。用双线性变换法，设计满足上述指标的数字ChebyshevI型低通滤波器。用一个图形窗口，包括三个子图，分析显示滤波器的幅频响应、以dB为纵坐标的幅频响应和相频响应。已知模拟信号其中，，，，取采样频率。以对进行取样，得到。用一个图形窗口，包括两个子图，分别显示以及()的波形；用FFT对进行谱分析，要求频率分辨率不超过5Hz。求出一个记录长度中的最少点数，并用一个图形窗口，包括两个子图，分别显示以及的幅值；用要求1中设计的线性相位低通数字滤波器对进行滤波，求出滤波器的输出，并用FFT对进行谱分析，要求频率分辨率不超过5Hz。求出一个记录长度中的最少点数，并用一个图形窗口，包括四个子图，分别显示（）、、（）和的幅值；用要求2中设计的Chebyshev低通数字滤波器对进行滤波，求出滤波器的输出，并用FFT对进行谱分析，要求频率分辨率不超过5Hz。求出一个记录长度中的最少点数，并用一个图形窗口，包括四个子图，分别显示（）、、（）和的幅值。三、设计思路1）窗函数设计数字线性相位FIR低通滤波器1、根据设计给出的技术指标，计算允许的过渡带宽度及阻带衰减，初步选择窗函数和N值。2、若选用理想低通逼近，则计算积分，求出理想低通的冲激响应hd(n);理想低通的截止频率选择为wc=(ws+wp)/2,hd(n)=sin[wc(n-a)]/pi*n-a);3、将hd(n)与窗函数相乘的FIR数字滤波器的冲激响应h(n)：h(n)=hd(n)w(n)4、计算FIR数字滤波器的频率响应，并验证是否达到所要求的指标。2）双线性变换法设计的数字ChebyshevI型低通滤波器1、根据设计给出的数字滤波器的技术指标，计算模拟原型低通的指标2、模拟切比雪夫1型滤波器的计算，利用function[b,a]=afd_cheb1(wp,ws,Rp,As)进行切比雪夫1型模拟滤波器的设计；利用function[b,a]=cheb1ap_o(N,Rp,omega)非归一化的切比雪夫1型模拟原型低通滤波器的设计3、双线性变换法将Ha(s)转换为H(z)4、直接型转换为级联型，计算数字滤波器的频率响应，并验证是否达到所要求的指标。3）模拟信号的波形产生和采样及FFT频谱分析和FIR滤波及IIR滤波N：一个记录长度中的取样数F：频率分辨率fs：采样周期最小记录长度的选择：由tp=NT=1/F可知，给定抽样频率和频率分辨率后，就可以得出最小记录长度N=fs/Ffft(x)：计算N点的DFT,N是序列的长度fft(x，L)：计算L点的DFT,若N>L，则将原序列截短为L点的序列，再计算其L点的DFT；若N<L，则将原序列补零为L点的序列，再计算其L点的DFT；ifft(x)：计算N点的IDFT,N是频域序列的长度ifft(x，L)：计算L点的IDFT,若N>L，则将原序列截短为L点的序列，再计算其L点的IDFT；若N<L，则将原序列补零为L点的序列，再计算其L点的IDFT；为了提高fft和ifft的计算效率，应尽量使序列长度N为2的幂次方，或将序列补零使其长度为2的幂次方。filter(b,a,x)：可采用此函数给出滤波器的输出b=H(z)的分子多项式系数（对FIRb=h)a=H(z)的分母多项式系数（对FIRa=[1])三、程序代码1、主程序代码%一、数字线性相位FIR低通滤波器%要求1

clear;

wp=0.26*pi;

ws=0.41*pi;

tr_width=ws-wp;%过度带宽

N=ceil(6.6*pi/tr_width)+1;%阶数

n=0:N-1;

wc=(ws+wp)/2;%理想低通的截止频率

hd=ideal_lp(wc,N);%理想低通滤波器的冲激响应

w_ham=(hamming(N))';%哈明窗的窗函数

h=hd.*w_ham;%FIR滤波器的冲激响应

[db,mag,pha,grd,w]=freqz_m(h,[1]);

delta_w=2*pi/1000;%取1000点

Rp=-(min(db(1:1:wp/delta_w+1)));%实际的通带衰减

As=-round(max(db(ws/delta_w+1:50)));%实际的最小阻带衰减

%作图

figure(1);

subplot(2,2,1);stem(n,hd);ylabel('hd(n)');xlabel('n');

title('单位冲激响应');grid;

axis([045-0.10.5]);

subplot(2,2,2);plot(w/pi,mag);title('幅度响应');grid;

axis([01.5-0.51.5]);xlabel('以\pi为单位的频率');ylabel('幅度');

subplot(2,2,3);plot(w/pi,db);title('幅度响应(dB)');grid;

axis([01.5-1600]);

xlabel('以为单位的频率');ylabel('对数幅度/dB');

subplot(2,2,4);plot(w/pi,pha);title('相位响应');grid;

axis([01-44]);

xlabel('以\pi为单位的频率');ylabel('相位');

%二、双线性变换法设计数字切比雪夫I型低通滤波器%要求2

%双线性变换法设计数字切比雪夫I型低通滤波器

%数字滤波器指标

wp=0.26*pi;ws=0.41*pi;Rp=0.75;As=50;

%转换为模拟原型低通指标

T=1;Fs=1/T;

omegap=(2/T)*tan(wp/2);

omegas=(2/T)*tan(ws/2);

%模拟切比雪夫I型滤波器得计算

[cs,ds]=afd_cheb1(omegap,omegas,Rp,As);

%双现性变换法

[b,a]=bilinear(cs,ds,Fs);

%计算数字滤波器的频率响应

[db,mag,pha,grd,w]=freqz_m(b,a);

figure(1);

subplot(3,1,1);plot(w/pi,mag);title('幅度响应');grid;

axis([00.801.1]);ylabel('幅度');xlabel('以\pi为单位的频率');

subplot(3,1,2);plot(w/pi,db);title('幅度响应(dB)');grid;

axis([00.8-601]);

xlabel('以\pi为单位的频率');ylabel('对数幅度/dB');

subplot(3,1,3);plot(w/pi,pha);title('相位响应');grid;

axis([00.8-44]);

ylabel('相位');xlabel('以\pi为单位的频率');

%三（1）、以对xt进行抽样，得到Xnn=0:1:511;fs=10000;t=0:0.000001:10;xt=2*sin(2*pi*120*t)+5*sin(2*pi*4980*t)+8*cos(2*pi*3250*t)+7.5*cos(2*pi*1150*t);T=1/fs;t0=n*T;%时间离散化Xn=2*sin(240*pi*t0)+5*sin(2*4980*pi*t0)+8*cos(2*pi*3250*t0)+7.5*cos(2*pi*1150*t0);figure(3)%显示以及()的波形；subplot(2,1,1);stem(n,Xn);title('Xn');grid;ylabel('Xn');xlabel('n');axis([01000-2525])subplot(2,1,2);plot(t,xt);title('xt曲线');grid;ylabel('xt');xlabel('t');axis([00.01-2525]);%三（2）、对Xn（0<=n<=Nx-1）进行FFT变换求Xk，对Xn进行谱分析%要求3--（2）

%信号的FFT分析

fs=10000;f(1)=120;f(2)=4980;f(3)=3250;f(4)=1150;

F=5;Nx=fs/F;T=1/fs;n=0:Nx-1;

Xn=2*sin(2*pi*f(1)*n/fs)+5*sin(2*pi*f(2)*n/fs)+8*cos(2*pi*f(3)*n/fs)+7.5*cos(2*pi*f(4)*n/fs);

fprintf('\n***直接FFT时一个记录长度最小点数Nx=%2.0f\n',N);

k=0:Nx-1;

Xk=fft(Xn);

magX=abs(Xk);%phaX=angle(Xk)figure(4)%显示以及的幅值；

subplot(2,1,1);stem(n,Xn);title('x(n)');grid;

ylabel('x(n)');xlabel('n');axis([0200-2525]);

subplot(2,1,2);stem(k,magX);title('x(k)');grid;

ylabel('magX');xlabel('k');axis([02048010000]);

%三（3）、用线性相位低通滤波器对Xn进行滤波得到y1n，对y1n进行FFT变换%要求3--(3)

Fy1=5;Ny1=fs/Fy1;

fprintf('\n**y1n用FFT进行谱分析时一个记录长度最小点数Ny1=%2.0f\n',Ny1);

y1n=filter(h,1,Xn);%Xn经过用哈明窗设计的线性相位低通滤波器滤波得到y1n

Y1k=fft(y1n);%对y1n进行FFT变换

magY1k=abs(Y1k);

figure(3);

subplot(2,2,1);stem(n,Xn);title('x(n)');grid;

ylabel('x(n)');xlabel('n');axis([0200-2525]);

subplot(2,2,2);stem(k,magX);title('X(k)');grid;

ylabel('magX');xlabel('k');axis([02048010000]);

subplot(2,2,3);stem(n,y1n);title('y1(n)');grid;

ylabel('y1(n)');xlabel('n');axis([0100-1010]);

subplot(2,2,4);stem(k,magY1k);title('Y1(k)');grid;

ylabel('magY1');xlabel('k');axis([0204808000]);

%三（4）、用由双线性变换法设计的数字切比雪夫I型低通滤波器对Xn进行滤波得到y2n，对y2n进行FFT变换%要求3--（4）

Ny2=fs/5;

fprintf('\n***y2n用FFT进行谱分析时一个记录长度最小点数Ny2=%2.0f\n',Ny2);

y2n=filter(b,a,Xn);%Xn经过用双线性变换法设计的数字切比雪夫I型低通滤波器滤波得到y1n

Yk2=fft(y2n);%对y2n进行FFT变换

magYk2=abs(Yk2);

figure(3);

subplot(2,2,1);stem(n,Xn);title('x(n)');grid;

ylabel('Xn');xlabel('n');axis([0200-2525]);

subplot(2,2,2);stem(k,magX);title('x(k)');grid;

ylabel('magX');xlabel('k');axis([02048010000]);

subplot(2,2,3);stem(n,y2n);title('y2(n)');grid;

ylabel('y2n');xlabel('n');axis([0100-1010]);

subplot(2,2,4);stem(k,magYk2);title('Y2(k)');grid;

ylabel('magY2');xlabel('k'