1、function xn = idft(Xk,N)WNnk=WN.(-nk);xn=(Xk*WNnk)/N;end%循环移位函数functiony=cirshift(x,m,N)if length(x)N error( N必须大于等于x(n)的长度)x=x zeros(1,N-length(x);n=0:N-1;n=mod(n-m,N);y=x(n+1);%计算序列x1(n) 和x2(n)之间的N点循环卷积函数functiony=circonvt(x1,x2,N)if length(x1) N必须大于等于x1的长度if length(x2) N必须大于等于x2的长度x1=x1 zeros(1,N
2、-length(x1);x2=x2 zeros(1,N-length(x2);m=0:x2=x2(mod(-m,N)+1);H=zeros(N,N);for n=1: H(n,:)=cirshift(x2,n-1,N);y=x1*H;例3.2(P139) x=1,1,1,1,1,1;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(
3、DTFT的幅度);xlabel(以pi为单位的频率ylabel(幅度 subplot(2,1,2);plot(w/pi,angX); axis(-1 1 -200 200);title(DTFT的相位 xlabel(相位 N=6; X=dft(x,N);phaX=angle(X)*180/pi; k=0: figure(2)stem(k,magX);DFT的幅度kstem(k,phaX);DFT的相位Xk = Columns 1 through 5 6.0000 -0.0000 - 0.0000i -0.0000 - 0.0000i 0 - 0.0000i 0.0000 - 0.0000i C
4、olumn 6 0.0000 - 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; magX=abs(X),angX=angle(X)stem(k,angX); 1.1611 1.2250 + 0.7289i 2.0522 + 3.5598i 0.9030 - 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 -
5、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:10;x=8*(0.8).n;N=11; y=cirshift(x,6,N);stem(n,x);序列x(n)nx(n)stem(n,y);x(n)的循环移位y(n)例3.5(P144) x1=1,2,2;x2=5,4,3,2,1;N=5; y=cir
6、convt(x1,x2,N)y =11 16 21 16 11例3.6(P144)x1=(0.8).n;6;x2=exp(-n);N=10; 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;r=n*T;x=r.*exp(r); X=fft(x,N); magX=abs(X),phaX=angle(X) subplot(3,1,1);axis(0 32 0 1); subplot(3,1,2);axis(0 32 0 8); subplot(3,1
7、,3);axis(0 32 -4 4); 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 through 30 0.2363 0.2475 0.2625 0.2824 0.3087 0.3441 0.3929 0.4630 0.5704 0.7518 Colu
8、mns 31 through 32 1.1182 2.2215phaX = 0 1.5747 1.7204 1.8347 1.9411 2.0443 2.1459 2.2467 2.3469 2.4467 2.5463 2.6457 2.7450 2.8442 2.9434 3.0425 3.1416 -3.0425 -2.9434 -2.8442 -2.7450 -2.6457 -2.5463 -2.4467 -2.3469 -2.2467 -2.1459 -2.0443 -1.9411 -1.8347 -1.7204 -1.5747例3.8(P148) N=21;L=256;f1=120;
9、f2=140;fs=400;ws=2*pi*fs;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);subplot(2,1,1);plot(w,abs(X);幅度谱频率(Hz) N=21axis(-200 200 0 15);subplot(2,1,2);axis(-200 200 0 25);频率(Hz) N=11例3.10(P152)20; x=cos(0.1*pi*n); h=(0.8).n; L=length(x)+length(h)-1; X=fft(x,L);
10、H=fft(h,L); y=ifft(X.*H)stem(n,h);序列h(n)L-1;stem(n,real(y); ylabel(卷积结果y(n) 1.0000 1.7511 2.2099 2.3557 2.1936 1.7548 1.0949 0.2881 -0.5785 -1.4139 -2.1311 -2.6559 -2.9338 -2.9348 -2.6569 -2.1255 -1.3914 -0.5253 0.3888 1.2621 2.0097 1.5985 1.2700 1.0086 0.8014 0.6383 0.5106 0.4114 0.3345 0.2751 Colu
11、mns 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.11(P153) b=0.0181,0.0543,0.0543,0.0181;a=1,-1.7600,1.1829,-0.2781;ch=impseq(0,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)
12、subplot(3,1,1);plot(n,x);axis(0 20 -1.5 1.5);信号X(n)subplot(3,1,2);axis(0 20 -0.5 0.5);系统冲激响应h(n)subplot(3,1,3);plot(n,y);x(n)滤波的结果y(n) 0.0130 0.0795 0.2283 0.4246 0.5909 0.6591 0.6050 0.4556 0.2503 -0.0104 -0.3476 -0.7126 -0.9934 -1.1208 -1.0877 -0.9072 -0.6222 -0.3122 -0.0386 0.2081 0.4825 0.7717 0
13、.9482 0.9129 0.7001 0.4154 0.1572 -0.0203 -0.1061 -0.1199 -0.0922 -0.0496 -0.0113 0.0121 0.0217 0.0204 0.0144 0.0063 0.0008 -0.0023 -0.0013四、实验分析根据实验绘出的图形,与实际运算出的结果相比较,可知,利用ifft函数与fft函数求出的值与实际求出的IDFT变换与DFT变换值的误差相差较小。五、实验结论如果采样数据过少,运用Fourier变换不能分辨其中的频率成分。添加零后可增加频谱的数据个数,谱的密度增高了,但仍不能分辨,其中的频率成分,即谱的分辨率没有
14、提高。只有数据点数足够多时才能分辨其中的频率成分。实验第4章matlab实现数字巴特沃斯滤波器设计一、实验目的1、掌握设计IIR滤波器的原理和具体设计方法;2、会在计算机上用IIR滤波器进行数字滤波;3、掌握对IIR滤波器特性分析。IIR数字滤波器相关知识。IIR数字滤波器的设计一般是利用目前已经很成熟的模拟滤波器的设计方法来进行设计,通常采用模拟滤波器原型有butterworth函数、chebyshev函数、bessel函数、椭圆滤波器函数等。%butterworth低通滤波器原型设计函数,要求WsWp0,AsRp0。function b,a=afd_butt(Wp,Ws,Rp,As)N=c
15、eil(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);%非归一化Butterworth模拟低通滤波器原形设计函数 %得到的b,a分别为传输函数分子、分母多项式系数;function b,a=u_buttap(N,Omegac);z,p,k=buttap(N); %归一化巴特沃思模拟低通滤
16、波器原形 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); %部分分式展开p=exp(p*T); %从模拟到数字极点对应关系 ,部分分式系数相同b,a=residuez(R,p,k); %将部分分式的形式变换成多项式之比的形式b=real(b %求出数字滤波器系数a=real(a%频率响应
17、函数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: %频率mag=abs(H); %响应幅度db=20*log10(mag+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; %以上步骤求出
18、系数 M=length(b); N=length(a);if NM b=b zeros(1,N-M);elseif M a=a zeros(1,M-N);else NM=0;K=floor(N/2); B=zeros(K,3); A=zeros(K,3);if K*2=N b=b 0; a=a 0;end broots=cplxpair(roots(b); %以下程序将每两个极点和两个零点组合成二阶因子aroots=cplxpair(roots(a); % roots:求多项式的根for i=1:2:2*K Brow=broots(i:i+1,: Brow=real(poly(Brow); B
19、(fix(i+1)/2,:)=Brow; Arow=aroots(i: Arow=real(poly(Arow); A(fix(i+1)/2,:)=Arow;%此函数产生理想低通滤波器的冲激响应function hd=ideal_lp(wc,M)alpha=(M-1)/2;(M-1);m=n-alpha+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_
20、invr(cs,ds,T)db,mag,pha,w=freqz_m(b,a);plot(w/pi,mag);digital filter Magnitude Responseaxis(0,1,0,1.1)plot(w/pi,db);digital filter Magnitude in DBaxis(0,1,-40,5);五、实验分析本次数字滤波器设计方法是基于MATLAB的数字滤波器的设计,是用学过的数字信号理论为依据,用MATLAB代码来实现,由滤波器的频谱图和滤波前后的信号的频谱图对比可知本设计选用双线性变换法设计的IIR滤波器比较好。在同样的技术指标的要求下,IIR滤波器所要求的阶数N也比较小,实现起来比较容易。六、实验结论在信号与信息的过滤、检测和预测等处理中,都要使用滤波器,数字滤波器是数字信号处理中使用最广泛的一种方法。IIR数字滤波器的设计过程中,可以借助模拟滤波器的设计成果或直接采用典型的滤波器类型,减少工作量.。通过这个实验,对设计数字滤波器的整个过程有了很好的掌握。
copyright@ 2008-2022 冰豆网网站版权所有
经营许可证编号:鄂ICP备2022015515号-1