五邑大学数字信号处理34章实验报告.docx

1、五邑大学数字信号处理34章实验报告五邑大学 数字信号处理课程设计报告题目：FFT和数字巴特沃斯滤波器的设计院 系 专 业 姓名学号 指导教师 实验第3章 matlab实现FFT应用一、实验目的熟悉matlab在数字信号处理中的应用，并掌握基二FFT算法的实现。二、实验原理FFT快速算法，matlabfft(x,N)实现应用，在MATLAB中，可以用函数X=fft（x，N）和x=ifft（X，N）计算N点序列的DFT正、反变换。N点序列的DFT和IDFT变换定义式如下： ，利用旋转因子具有周期性，可以得到快速算法（FFT）。三、程序代码与实验波形%计算N点的DFTfunction Xk = df

2、t(xn,N)n=0:1:N-1;k=0:1:N-1;WN=exp(-j*2*pi/N);nk=n*k;WNnk=WN.nk;Xk=xn*WNnkEnd%计算N点的逆DFTfunction xn = idft(Xk,N)n=0:1:N-1;k=0:1:N-1;WN=exp(-j*2*pi/N);nk=n*k;WNnk=WN.(-nk);xn=(Xk*WNnk)/N;end%循环移位函数functiony=cirshift(x,m,N)if length(x)N error( N必须大于等于x(n)的长度)endx=x zeros(1,N-length(x);n=0:N-1;n=mod(n-m,

3、N);y=x(n+1);%计算序列x1(n) 和x2(n)之间的N点循环卷积函数functiony=circonvt(x1,x2,N)if length(x1)N error( N必须大于等于x1的长度)endif length(x2)N error( N必须大于等于x2的长度)endx1=x1 zeros(1,N-length(x1);x2=x2 zeros(1,N-length(x2);m=0:1:N-1;x2=x2(mod(-m,N)+1);H=zeros(N,N);for n=1:1:N H(n,:)=cirshift(x2,n-1,N);endy=x1*H;例3.2（P139） x=

4、1,1,1,1,1,1;n=0:5; k=-200:200;w=(pi/100)*k; X=x*(exp(-j*pi/100).(n*k); magX=abs(X);angX=angle(X)*180/pi; figure(1) subplot(2,1,1);plot(w/pi,magX); axis(-1 1 0 6);grid; title(DTFT的幅度);xlabel(以pi为单位的频率);ylabel(幅度); subplot(2,1,2);plot(w/pi,angX);grid; axis(-1 1 -200 200);title(DTFT的相位); xlabel(以pi为单位的

5、频率);ylabel(相位); N=6; X=dft(x,N); magX=abs(X);phaX=angle(X)*180/pi; k=0:5; figure(2) subplot(2,1,1);stem(k,magX); title(DFT的幅度);xlabel(k); subplot(2,1,2);stem(k,phaX); title(DFT的相位);xlabel(k);Xk = Columns 1 through 5 6.0000 -0.0000 - 0.0000i -0.0000 - 0.0000i 0 - 0.0000i 0.0000 - 0.0000i Column 6 0.0

6、000 - 0.0000i例3.3（P141） n=0:9;x=0.8*cos(0.47*pi*n)+0.4*cos(0.53*pi*n); N=10; X=dft(x,N); magX=abs(X),angX=angle(X) k=0:9; subplot(2,1,1);stem(k,magX);title(DFT的幅度); xlabel(k); subplot(2,1,2);stem(k,angX);title(DFT的相位); xlabel(k);Xk = Columns 1 through 5 1.1611 1.2250 + 0.7289i 2.0522 + 3.5598i 0.903

7、0 - 2.9037i 0.8213 - 0.6744i Columns 6 through 10 0.8360 - 0.0000i 0.8213 + 0.6744i 0.9030 + 2.9037i 2.0522 - 3.5598i 1.2250 - 0.7289imagX = 1.1611 1.4254 4.1090 3.0409 1.0627 0.8360 1.0627 3.0409 4.1090 1.4254angX = 0 0.5367 1.0478 -1.2693 -0.6875 -0.0000 0.6875 1.2693 -1.0478 -0.5367例3.4（P143） n=0

8、:10;x=8*(0.8).n;N=11; y=cirshift(x,6,N); subplot(2,1,1);stem(n,x);title(序列x(n); xlabel(n);ylabel(x(n); subplot(2,1,2);stem(n,y);title(x(n)的循环移位); xlabel(n);ylabel(y(n);例3.5（P144） x1=1,2,2;x2=5,4,3,2,1;N=5; y=circonvt(x1,x2,N)y =11 16 21 16 11例3.6（P144） n=0:9;x1=(0.8).n; n=0:6;x2=exp(-n);N=10; y=circ

9、onvt(x1,x2,N)y = 1.0905 1.2008 1.0815 0.9096 0.7440 0.6012 0.4832 0.3866 0.3093 0.2474例3.7（P146） N=32;fs=100;T=1/fs; n=0:N-1;r=n*T;x=r.*exp(r); X=fft(x,N); magX=abs(X),phaX=angle(X) subplot(3,1,1);stem(n,x);axis(0 32 0 1); xlabel(n);title(序列x(n); k=0:N-1; subplot(3,1,2);stem(k,magX);axis(0 32 0 8);

10、xlabel(k);ylabel(DFT的幅度); subplot(3,1,3);stem(k,phaX);axis(0 32 -4 4); xlabel(k);ylabel(DFT的相位);magX = Columns 1 through 10 6.1357 2.2215 1.1182 0.7518 0.5704 0.4630 0.3929 0.3441 0.3087 0.2824 Columns 11 through 20 0.2625 0.2475 0.2363 0.2281 0.2226 0.2194 0.2183 0.2194 0.2226 0.2281 Columns 21 thr

11、ough 30 0.2363 0.2475 0.2625 0.2824 0.3087 0.3441 0.3929 0.4630 0.5704 0.7518 Columns 31 through 32 1.1182 2.2215phaX = Columns 1 through 10 0 1.5747 1.7204 1.8347 1.9411 2.0443 2.1459 2.2467 2.3469 2.4467 Columns 11 through 20 2.5463 2.6457 2.7450 2.8442 2.9434 3.0425 3.1416 -3.0425 -2.9434 -2.8442

12、 Columns 21 through 30 -2.7450 -2.6457 -2.5463 -2.4467 -2.3469 -2.2467 -2.1459 -2.0443 -1.9411 -1.8347 Columns 31 through 32 -1.7204 -1.5747例3.8（P148） N=21;L=256;f1=120;f2=140;fs=400;T=1/fs;ws=2*pi*fs;n=0:N-1;x=cos(2*pi*f1*n*T)+cos(2*pi*f2*n*T);X=fftshift(fft(x,L);w=(-ws/2+(0:L-1)*ws/L)/(2*pi);subpl

13、ot(2,1,1);plot(w,abs(X);ylabel(幅度谱);xlabel(频率(Hz) N=21);axis(-200 200 0 15);N=11;n=0:N-1;x=cos(2*pi*f1*n*T)+cos(2*pi*f2*n*T);X=fftshift(fft(x,L);subplot(2,1,2);plot(w,abs(X);axis(-200 200 0 25);ylabel(幅度谱);xlabel(频率(Hz) N=11);例3.10（P152） n=0:20; x=cos(0.1*pi*n); h=(0.8).n; L=length(x)+length(h)-1; X

14、=fft(x,L); H=fft(h,L); y=ifft(X.*H) subplot(3,1,1);stem(n,x);ylabel(序列x(n);xlabel(n); subplot(3,1,2);stem(n,h);ylabel(序列h(n);xlabel(n); n=0:L-1; subplot(3,1,3);stem(n,real(y); ylabel(卷积结果y(n);xlabel(n);y = Columns 1 through 10 1.0000 1.7511 2.2099 2.3557 2.1936 1.7548 1.0949 0.2881 -0.5785 -1.4139 C

15、olumns 11 through 20 -2.1311 -2.6559 -2.9338 -2.9348 -2.6569 -2.1255 -1.3914 -0.5253 0.3888 1.2621 Columns 21 through 30 2.0097 1.5985 1.2700 1.0086 0.8014 0.6383 0.5106 0.4114 0.3345 0.2751 Columns 31 through 40 0.2288 0.1923 0.1626 0.1375 0.1155 0.0952 0.0762 0.0581 0.0410 0.0254 Column 410.0115例3

16、.11（P153） b=0.0181,0.0543,0.0543,0.0181;a=1,-1.7600,1.1829,-0.2781;ch=impseq(0,0,20);n=0:20;h=filter(b,a,ch);x=cos(0.1*pi*n)+0.32*randn(size(n);L=length(x)+length(h)-1;X=fft(x,L);H=fft(h,L);y=ifft(X.*H)subplot(3,1,1);plot(n,x);axis(0 20 -1.5 1.5);ylabel(信号X(n);xlabel(n);subplot(3,1,2);stem(n,h);axis

17、(0 20 -0.5 0.5);ylabel(系统冲激响应h(n);xlabel(n);n=0:L-1;subplot(3,1,3);plot(n,y);axis(0 20 -1.5 1.5);ylabel(x(n)滤波的结果y(n);xlabel(n);y = Columns 1 through 10 0.0130 0.0795 0.2283 0.4246 0.5909 0.6591 0.6050 0.4556 0.2503 -0.0104 Columns 11 through 20 -0.3476 -0.7126 -0.9934 -1.1208 -1.0877 -0.9072 -0.622

18、2 -0.3122 -0.0386 0.2081 Columns 21 through 30 0.4825 0.7717 0.9482 0.9129 0.7001 0.4154 0.1572 -0.0203 -0.1061 -0.1199 Columns 31 through 40 -0.0922 -0.0496 -0.0113 0.0121 0.0217 0.0204 0.0144 0.0063 0.0008 -0.0023 Column 41 -0.0013四、实验分析根据实验绘出的图形，与实际运算出的结果相比较，可知，利用ifft函数与fft函数求出的值与实际求出的IDFT变换与DFT变

19、换值的误差相差较小。五、实验结论如果采样数据过少，运用Fourier变换不能分辨其中的频率成分。添加零后可增加频谱的数据个数，谱的密度增高了，但仍不能分辨，其中的频率成分，即谱的分辨率没有提高。只有数据点数足够多时才能分辨其中的频率成分。实验第4章matlab实现数字巴特沃斯滤波器设计一、实验目的1、掌握设计IIR滤波器的原理和具体设计方法；2、会在计算机上用IIR滤波器进行数字滤波；3、掌握对IIR滤波器特性分析。二、实验原理IIR数字滤波器相关知识。IIR数字滤波器的设计一般是利用目前已经很成熟的模拟滤波器的设计方法来进行设计，通常采用模拟滤波器原型有butterworth函数、cheby

20、shev函数、bessel函数、椭圆滤波器函数等。三、程序代码与实验波形%butterworth低通滤波器原型设计函数，要求WsWp0，AsRp0。function b,a=afd_butt(Wp,Ws,Rp,As)N=ceil(log10(10(Rp/10)-1)/(10(As/10)-1)/(2*log10(Wp/Ws);%ceil 朝正无穷大方向取整；fprintf(n Butterworth Filter Order=%2.0fn,N)OmegaC=Wp/(10(Rp/10)-1)(1/(2*N) %求对应于N的3db截止频率；b,a=u_buttap(N,OmegaC);%非归一化B

21、utterworth模拟低通滤波器原形设计函数 %得到的b,a分别为传输函数分子、分母多项式系数；function b,a=u_buttap(N,Omegac);z,p,k=buttap(N); %归一化巴特沃思模拟低通滤波器原形 p=p*Omegac; %将 代入上式，相当于分子乘以 ，极点乘以 k=k*OmegacN;B=real(poly(z); %poly为构造具有指定根的多项式 real为求实部b=k*B;a=real(poly(p); %利用脉冲响应不变法从模拟到数字滤波器变换函数function b,a=imp_invr(c,d,T)R,p,k=residue(c,d); %部分

22、分式展开p=exp(p*T); %从模拟到数字极点对应关系 ,部分分式系数相同b,a=residuez(R,p,k); %将部分分式的形式变换成多项式之比的形式b=real(b); %求出数字滤波器系数a=real(a);%频率响应函数freqz的修正，此函数可获得滤波器的幅值响应、相位响应及群延迟响应function db,mag,pha,w=freqz_m(b,a)H,w=freqz(b,a,1000,whole); %在0-2*pi之间选取N个点计算频率响应H=(H(1:501); %频率响应 w=(w(1:501); %频率mag=abs(H); %响应幅度db=20*log10(ma

23、g+eps)/max(mag); %增益pha=angle(H); %相位 %变直接形式为级联形式function b0,B,A=dir2cas(b,a)b0=b(1);b=b/b0;a0=a(1);a=a/a0;b0=b0/a0; %以上步骤求出系数 M=length(b); N=length(a);if NM b=b zeros(1,N-M);elseif MN a=a zeros(1,M-N);else NM=0;endK=floor(N/2); B=zeros(K,3); A=zeros(K,3);if K*2=N b=b 0; a=a 0;end broots=cplxpair(ro

24、ots(b); %以下程序将每两个极点和两个零点组合成二阶因子aroots=cplxpair(roots(a); % roots：求多项式的根for i=1:2:2*K Brow=broots(i:1:i+1,:); Brow=real(poly(Brow); B(fix(i+1)/2,:)=Brow; Arow=aroots(i:1:i+1,:); Arow=real(poly(Arow); A(fix(i+1)/2,:)=Arow;end%此函数产生理想低通滤波器的冲激响应function hd=ideal_lp(wc,M)alpha=(M-1)/2; n=0:(M-1); m=n-alp

25、ha+eps; hd=sin(wc*m)./(pi*m); %利用脉冲响应不变法设计巴特沃思滤波器wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;T=1;OmegaP=wp/T;OmegaS=ws/T;cs,ds=afd_butt(OmegaP,OmegaS,Rp,As);b,a=imp_invr(cs,ds,T)db,mag,pha,w=freqz_m(b,a);subplot(2,1,1);plot(w/pi,mag);title(digital filter Magnitude Response)axis(0,1,0,1.1)subplot(2,1,2);plot(w/pi

26、,db);title(digital filter Magnitude in DB)axis(0,1,-40,5);五、实验分析本次数字滤波器设计方法是基于MATLAB的数字滤波器的设计，是用学过的数字信号理论为依据，用MATLAB代码来实现，由滤波器的频谱图和滤波前后的信号的频谱图对比可知本设计选用双线性变换法设计的IIR滤波器比较好。在同样的技术指标的要求下，IIR滤波器所要求的阶数N也比较小，实现起来比较容易。六、实验结论在信号与信息的过滤、检测和预测等处理中，都要使用滤波器，数字滤波器是数字信号处理中使用最广泛的一种方法。IIR数字滤波器的设计过程中，可以借助模拟滤波器的设计成果或直接采用典型的滤波器类型，减少工作量.。通过这个实验，对设计数字滤波器的整个过程有了很好的掌握。