2017年西电电院数字信号处理上机实验报告五
西电DSP实验报告

西安电子科技大学数字信号处理(DSP)课程实验报告实验名称 VISUAL DSP++的使用入门电子工程学院 1402071 班同作者崔健孟智超杨伟祺实验日期 2017 年 5 月 16 日实验一:VISUAL DSP++的使用入门一、实验目的:熟悉VISUAL DSP++的开发环境二、实验内容:练习一:启动Visual DSP++,建立一个用C源代码的工程(Project),同时用调试器来评估用C语言所编写代码的性能;练习二:创立一个新的工程,修改源码来调用一个汇编(asm)程序,重新编译工程,用调试器来评估用汇编语言所写程序的性能;练习三:利用调试器的绘图(plot)功能来图形显示一个卷积算法中的多个数据的波形;练习四:利用调试器的性能统计功能(Statistical profile)来检查练习三中卷积算法的效率。
利用所收集到的性能统计数据就能看出算法中最耗时的地方。
三、实验步骤及实验结果:练习一:1)新建工程进入 Visual DSP++,显示Visual DSP++的集成开发和调试环境窗口,选择菜单File 中Open 打开文件:…\unit_1\dot_product_c \dotprodc.dpj。
2)编译 dotprodc工程在菜单 Project中选择 Build Project来对工程进行编译。
在本例子中,编译器会检测到一个未定义的错误,显示为:“.\dotprod_main.c”,line 115:error #20:identifier“itn”is undefined itn i;将该错误改正后,保存并重新编译,没有错误出现,编译成功。
3)运行VsualDSP++调试器在编译完成后,环境将自动进入调试状态,对于初次进入debugger,将显示对象选择对话框,在其中指定对象和处理器信息。
4)运行dotprod.c从 Debug菜单中选择 Run项,程序将被执行,其输出结果在 Output window中显示。
西安电子科技大学数字信号处理上机作业

数字信号处理MATLAB上机作业M 2.21.题目The square wave and the sawtooth wave are two periodic sequences as sketched in figure ing the function stem. The input data specified by the user are: desired length L of the sequence, peak value A, and the period N. For the square wave sequence an additional user-specified parameter is the duty cycle, which is the percent of the period for which the signal is positive. Using this program generate the first 100 samples of each of the above sequences with a sampling rate of 20 kHz ,a peak value of 7, a period of 13 ,and a duty cycle of 60% for the square wave.2.程序% 用户定义各项参数参数A = input('The peak value =');L = input('Length of sequence =');N = input('The period of sequence =');FT = input('The desired sampling frequency =');DC = input('The square wave duty cycle = ');% 产生所需要的信号t = 0:L-1;T = 1/FT;x = A*sawtooth(2*pi*t/N);y = A*square(2*pi*(t/N),DC);% Plotsubplot(2,1,1)stem(t,x);ylabel('幅度');xlabel('n');subplot(2,1,2)stem(t,y);ylabel('幅度');xlabel('n');3.结果4.结果分析M 2.41.题目(a)Write a matlab program to generate a sinusoidal sequence x[n]= Acos(ω0 n+Ф) and plot thesequence using the stem function. The input data specified by the user are the desired length L, amplitude A, the angular frequency ω0 , and the phase Фwhere 0<ω0 <pi and 0<=Ф<=2pi. Using this program generate the sinusoidal sequences shown in figure 2.15. (b)Generate sinusoidal sequences with the angular frequencies given in Problem 2.22.Determine the period of each sequence from the plot and verify the result theoretically. 2.程序%用户定义的参数L = input('Desired length = ');A = input('Amplitude = ');omega = input('Angular frequency = ');phi = input('Phase = ');%信号产生n = 0:L-1;x = A*cos(omega*n + phi);stem(n,x);xlabel('n');ylabel('幅度');title(['\omega_{o} = ',num2str(omega)]);3.结果(a)ω0=0ω0=0.1πω0=0.8πω0=1.2π(b)ω0=0.14πω0=0.24πω0=0.34πω0=0.68πω0=0.75π4.结果分析M 2.51.题目Generate the sequences of problem 2.21(b) to 2.21(e) using matlab.2.程序(b)n = 0 : 99;x=sin(0.6*pi*n+0.6*pi);stem(n,x);xlabel('n');ylabel('幅度');(c)n = 0 : 99;x=2*cos(1.1*pi*n-0.5*pi)+2*sin(0.7*pi*n);stem(n,x);xlabel('n');ylabel('幅度');(d)n = 0 : 99;x=3*sin(1.3*pi*n-4*cos(0.3*pi*n+0.45*pi));stem(n,x);xlabel('n');ylabel('幅度');(e)n = 0 : 99;x=5*sin(1.2*pi*n+0.65*pi)+4*sin(0.8*pi*n)-cos(0.8*pi*n);stem(n,x);xlabel('n');ylabel('幅度');(f)n = 0 : 99;x=mod(n,6);stem(n,x);xlabel('n');ylabel('幅度');3.结果(b)(c)(d)(e)(f)4.结果分析M 2.61.题目Write a matlab program to plot a continuous-time sinusoidal signal and its sampled version and verify figure 2.19. You need to use the hold function to keep both plots.2.程序%用户定义的参数fo = input('Frequency of sinusoid in Hz = ');FT = input('Samplig frequency in Hz = ');%产生信号t = 0:0.001:1;g1 = cos(2*pi*fo*t);plot(t,g1,'-')xlabel('时间t');ylabel('幅度')holdn = 0:1:FT;gs = cos(2*pi*fo*n/FT);plot(n/FT,gs,'o');hold off3.结果4.结果分析M 3.11.题目Using program 3_1 determine and plot the real and imaginary parts and the magnitude and phase spectra of the following DTFT for various values of r and θ:G(e jω)=1, 0<r<1.1−2r(cosθ)e−jω+r2e−2jω2.程序%program 3_1%discrete-time fourier transform computatition%k=input('Number of frequency points = ');num=input('Numerator coefficients= ');den=input('Denominator coefficients= ');%computer the frequency responsew=0:pi/k:pi;h=freqz(num,den,w);%plot the frequency responsesubplot(221)plot(w/pi,real(h));gridtitle('real part')xlabel('\omega/\pi');ylabel('Amplitude') subplot(222)plot(w/pi,imag(h));gridtitle('imaginary part')xlabel('\omega/\pi');ylabel('Amplitude') subplot(223)plot(w/pi,abs(h));gridtitle('magnitude spectrum')xlabel('\omega/\pi');ylabel('magnitude') subplot(224)plot(w/pi,angle(h));gridtitle('phase spectrum')xlabel('\omega/\pi');ylabel('phase,radians')3.结果(a)r=0.8 θ=π/6(b)r=0.6 θ=π/34.结果分析M 3.41.题目Using matlab verify the following general properties of the DTFT as listed in Table 3.2:(a) Linearity, (b) time-shifting, (c) frequency-shifting, (d) differentiation-in-frequency, (e) convolution, (f) modulation, and (g) Parseval’s relation. Since all data in matlab have to be finite-length vectors, the sequences to be used to verify the properties are thus restricted to be of finite length.2.程序%先定义两个信号N = input('The length of the sequence = ');k = 0:N-1;%g为正弦信号g = 2*sin(2*pi*k/(N/2));%h为余弦信号h = 3*cos(2*pi*k/(N/2));[G,w] = freqz(g,1);[H,w] = freqz(h,1);%*************************************************************************%% 线性性质alpha = 0.5;beta = 0.25;y = alpha*g+beta*h;[Y,w] = freqz(y,1);figure(1);subplot(211),plot(w/pi,abs(Y));xlabel('\omega/\pi');ylabel('|Y(e^j^\omega)|');title('线性叠加后的频率特性');grid;% 画出Y 的频率特性subplot(212),plot(w/pi,alpha*abs(G)+beta*abs(H));xlabel('\omega/\pi');ylabel('\alpha|G(e^j^\omega)|+\beta|H(e^j^\omega)|');title('线性叠加前的频率特性');grid;% 画出alpha*G+beta*H 的频率特性%*************************************************************************% % 时移性质n0 = 10;%时移10个的单位y2 = [zeros([1,n0]) g];[Y2,w] = freqz(y2,1);G0 = exp(-j*w*n0).*G;figure(2);subplot(211),plot(w/pi,abs(G0));xlabel('\omega/\pi');ylabel('|G0(e^j^\omega)|');title('G0的频率特性');grid;% 画出G0的频率特性subplot(212),plot(w/pi,abs(Y2));xlabel('\omega/\pi');ylabel('|Y2(e^j^\omega)|');title('Y2的频率特性');grid;% 画出Y2 的频率特性%*************************************************************************% % 频移特性w0 = pi/2; % 频移pi/2r=256; %the value of w0 in terms of number of samplesk = 0:N-1;y3 = g.*exp(j*w0*k);[Y3,w] = freqz(y3,1);% 对采样的512个点分别进行减少pi/2,从而生成G(exp(w-w0))k = 0:511;w = -w0+pi*k/512;G1 = freqz(g,1,w);figure(3);subplot(211),plot(w/pi,abs(Y3));xlabel('\omega/\pi');ylabel('|Y3(e^j^\omega)|');title('Y3的频率特性');grid;% 画出Y3的频率特性subplot(212),plot(w/pi,abs(G1));xlabel('\omega/\pi');ylabel('|G1(e^j^\omega)|');title('G1的频率特性');grid;% 画出G1 的频率特性%*************************************************************************% % 频域微分k = 0:N-1;y4 = k.*g;[Y4,w] = freqz(y4,1);%在频域进行微分y0 = ((-1).^k).*g;G2 = [G(2:512)' sum(y0)]';delG = (G2-G)*512/pi;figure(4);subplot(211),plot(w/pi,abs(Y4));xlabel('\omega/\pi');ylabel('|Y4(e^j^\omega)|');title('Y4的频率特性');grid;% 画出Y4的频率特性subplot(212),plot(w/pi,abs(delG));xlabel('\omega/\pi');ylabel('|delG(e^j^\omega)|');title('delG的频率特性');grid;% 画出delG的频率特性%*************************************************************************% % 相乘性质y5 = conv(g,h);%时域卷积[Y5,w] = freqz(y5,1);figure(5);subplot(211),plot(w/pi,abs(Y5));xlabel('\omega/\pi');ylabel('|Y5(e^j^\omega)|');title('Y5的频率特性');grid;% 画出Y5的频率特性subplot(212),plot(w/pi,abs(G.*H));%频域乘积xlabel('\omega/\pi');ylabel('|G.*H(e^j^\omega)|');title('G.*H的频率特性');grid;% 画出G.*H的频率特性%*************************************************************************% % 帕斯瓦尔定理y6 = g.*h;%对于freqz函数,在0到2pi直接取样[Y6,w] = freqz(y6,1,512,'whole');[G0,w] = freqz(g,1,512,'whole');[H0,w] = freqz(h,1,512,'whole');% Evaluate the sample value at w = pi/2% and verify with Y6 at pi/2H1 = [fliplr(H0(1:129)') fliplr(H0(130:512)')]';val = 1/(512)*sum(G0.*H1);% Compare val with Y6(129) i.e sample at pi/2 % Can extend this to other points similarly% Parsevals theoremval1 = sum(g.*conj(h));val2 = sum(G0.*conj(H0))/512;% Comapre val1 with val23.结果(a)(b)(c)(d)(e)4.结果分析M 3.81.题目Using matlab compute the N-point DFTs of the length-N sequences of Problem 3.12 for N=3, 5, 7, and 10. Compare your results with that obtained by evaluating the DTFTs computed in Problem 3.12 at ω= 2pik/N, k=0, 1,……N-1.2.程序%用户定义N的长度N = input('The value of N = ');k = -N:N;y1 = ones([1,2*N+1]);w = 0:2*pi/255:2*pi;Y1 = freqz(y1, 1, w);%对y1做傅里叶变换Y1dft = fft(y1);k = 0:1:2*N;plot(w/pi,abs(Y1),k*2/(2*N+1),abs(Y1dft),'o');grid;xlabel('归一化频率');ylabel('幅度');(a)clf;N = input('The value of N = ');k = -N:N;y1 = ones([1,2*N+1]);w = 0:2*pi/255:2*pi;Y1 = freqz(y1, 1, w);Y1dft = fft(y1);k = 0:1:2*N;plot(w/pi,abs(Y1),k*2/(2*N+1),abs(Y1dft),'o');xlabel('Normalized frequency');ylabel('Amplitude');(b)%用户定义N的长度N = input('The value of N = ');k = -N:N;y1 = ones([1,2*N+1]);y2 = y1 - abs(k)/N;w = 0:2*pi/255:2*pi;Y2 = freqz(y2, 1, w);%对y1做傅里叶变换Y2dft = fft(y2);k = 0:1:2*N;plot(w/pi,abs(Y2),k*2/(2*N+1),abs(Y2dft),'o');grid;xlabel('归一化频率');ylabel('幅度');(c)%用户定义N的长度N = input('The value of N = ');k = -N:N;y3 =cos(pi*k/(2*N));w = 0:2*pi/255:2*pi;Y3 = freqz(y3, 1, w);%对y1做傅里叶变换Y3dft = fft(y3);k = 0:1:2*N;plot(w/pi,abs(Y3),k*2/(2*N+1),abs(Y3dft),'o');grid;xlabel('归一化频率');ylabel('幅度');3.结果(a)N=3N=5 N=7N=10 (b)N=3N=5 N=7N=10 (c)N=3N=5 N=7N=104.结果分析M 3.191.题目Using Program 3_10 determine the z-transform as a ratio of two polynomials in z-1 from each of the partial-fraction expansions listed below:(a)X1(z)=−2+104+z−1−82+z−1,|z|>0.5,(b)X2(z)=3.5−21−0.5z−1−3+z−11−0.25z−2,|z|>0.5,(c)X3(z)=5(3+2z−1)2−43+2z−1+31+0.81z−2,|z|>0.9,(d)X4(z)=4+105+2z−1+z−16+5z−1+z−2,|z|>0.5.2.程序% Program 3_10% Partical-Fraction Expansion to rational z-Transform %r = input('Type in the residues = ');p = input('Type in the poles = ');k = input('Type in the constants = ');[num, den] = residuez(r,p,k);disp('Numberator polynominal coefficients');disp(num) disp('Denominator polynomial coefficients'); disp(den)4.结果分析M 4.61.题目Plot the magnitude and phase responses of the causal IIR digital transfer functionH(z)=0.0534(1+z−1)(1−1.0166z−1+z−2) (1−0.683z−1)(1−1.4461z−1+0.7957z−2).What type of filter does this transfer function represent? Determine the difference equation representation of the above transfer function.2.程序b=[0.0534 -0.00088644 -0.00088644 0.0534];a=[1 -2.1291 1.7833863 -0.5434631];figure(1)freqz(b,a);figure(2)[H,w]=freqz(b,a);plot(w/pi,abs(H)),grid;xlabel('Normalized Frequency (\times\pi rad/sample)'),ylabel('Magnitude');幅度化成真值之后:4.结果分析H(z)=0.0534−0.00088644z−1−0.00088644z−2+0.0534z−31−2.1291z−1+1.7833863z−2−0.5434631z−3M 4.71.题目Plot the magnitude and phase responses of the causal IIR digital transfer functionH(z)=(1−z−1)4(1−1.499z−1+0.8482z−2)(1−1.5548z−1+0.6493z−2).2.程序b=[1 -4 6 -4 1];a=[1 -3.0538 3.8227 -2.2837 0.5472]; figure(1)freqz(b,a);figure(2)[H,w]=freqz(b,a);plot(w/pi,abs(H)),grid;xlabel('Normalized Frequency (\times\pi rad/sample)'), ylabel('Magnitude');3.结果4.结果分析。
(完整word版)数字信号处理上机实验答案(第三版,第十章)

第十章 上机实验数字信号处理是一门理论和实际密切结合的课程,为深入掌握课程内容,最好在学习理论的同时,做习题和上机实验。
上机实验不仅可以帮助读者深入的理解和消化基本理论,而且能锻炼初学者的独立解决问题的能力。
本章在第二版的基础上编写了六个实验,前五个实验属基础理论实验,第六个属应用综合实验。
实验一 系统响应及系统稳定性。
实验二 时域采样与频域采样。
实验三 用FFT 对信号作频谱分析。
实验四 IIR 数字滤波器设计及软件实现。
实验五 FIR 数字滤波器设计与软件实现实验六 应用实验——数字信号处理在双音多频拨号系统中的应用任课教师根据教学进度,安排学生上机进行实验。
建议自学的读者在学习完第一章后作实验一;在学习完第三、四章后作实验二和实验三;实验四IIR 数字滤波器设计及软件实现在。
学习完第六章进行;实验五在学习完第七章后进行。
实验六综合实验在学习完第七章或者再后些进行;实验六为综合实验,在学习完本课程后再进行。
10.1 实验一: 系统响应及系统稳定性1.实验目的(1)掌握 求系统响应的方法。
(2)掌握时域离散系统的时域特性。
(3)分析、观察及检验系统的稳定性。
2.实验原理与方法在时域中,描写系统特性的方法是差分方程和单位脉冲响应,在频域可以用系统函数描述系统特性。
已知输入信号可以由差分方程、单位脉冲响应或系统函数求出系统对于该输入信号的响应,本实验仅在时域求解。
在计算机上适合用递推法求差分方程的解,最简单的方法是采用MA TLAB 语言的工具箱函数filter 函数。
也可以用MATLAB 语言的工具箱函数conv 函数计算输入信号和系统的单位脉冲响应的线性卷积,求出系统的响应。
系统的时域特性指的是系统的线性时不变性质、因果性和稳定性。
重点分析实验系统的稳定性,包括观察系统的暂态响应和稳定响应。
系统的稳定性是指对任意有界的输入信号,系统都能得到有界的系统响应。
或者系统的单位脉冲响应满足绝对可和的条件。
数字信号处理上机实验报告

实验一系统响应及系统稳定性一、实验目的(1)掌握求系统响应的方法。
(2)掌握时域离散系统的时域特性。
(3)分析、观察及检验系统的稳定性。
二、实验内容(1)给定一个低通滤波器的差分方程为y(n)=(n)+(n-1)+(n-1)输入信号x1(n)=R8(n)x2(n)=u(n)(a) 分别求出系统对x1(n)=R8(n) 和x2(n)=u(n)的响应序列,并画出其波形。
(b) 求出系统的单位冲响应,画出其波形。
实验程序:A=[1,];B=[,]; %%系统差分方程系数向量 B 和 Ax1n=[1 1 1 1 1 1 1 1 zeros(1,50)]; %产生信号 x1(n)=R8(n) x2n=ones(1,128); %产生信号 x2(n)=u(n)y1n=filter(B,A,x1n); %求系统对 x1(n)的响应 y1(n)n=0:length(y1n)-1;subplot(2,2,1);stem(n,y1n,'.');title('(a) 系统对 R_8(n)的响应y_1(n)');xlabel('n');ylabel('y_1(n)');y2n=filter(B,A,x2n); %求系统对 x2(n)的响应 y2(n)n=0:length(y2n)-1;subplot(2,2,2);stem(n,y2n,'.');title('(b) 系统对 u(n)的响应y_2(n)');xlabel('n');ylabel('y_2(n)');hn=impz(B,A,58); %求系统单位脉冲响应 h(n) n=0:length(hn)-1;subplot(2,2,3);y=hn;stem(n,hn,'.');title('(c) 系统单位脉冲响应h(n)');xlabel('n');ylabel('h(n)');运行结果图:(2)给定系统的单位脉冲响应为h1(n)=R10(n)h2(n)= δ(n)+δ(n-1)+δ(n-2)+δ(n-3)用线性卷积法分别求系统h1(n)和h2(n)对x1(n)=R8(n)的输出响应,并画出波形。
西安电子科技大学数字信号处理上机报告

数字信号处理大作业院系:电子工程学院学号:02115043姓名:邱道森实验一:信号、系统及系统响应一、实验目的(1) 熟悉连续信号经理想采样前后的频谱变化关系,加深对时域采样定理的理解。
(2) 熟悉时域离散系统的时域特性。
(3) 利用卷积方法观察分析系统的时域特性。
(4) 掌握序列傅里叶变换的计算机实现方法,利用序列的傅里叶变换对连续信号、离散信号及系统响应进行频域分析。
二、实验原理采样是连续信号数字处理的第一个关键环节。
对连续信号()a x t 进行理想采样的过程可用(1.1)式表示:()()()ˆa a xt x t p t =⋅ 其中()t xa ˆ为()a x t 的理想采样,()p t 为周期冲激脉冲,即 ()()n p t t nT δ∞=-∞=-∑()t xa ˆ的傅里叶变换()j a X Ω为 ()()s 1ˆj j j a a m X ΩX ΩkΩT ∞=-∞=-∑进行傅里叶变换,()()()j ˆj e d Ωt a a n X Ωx t t nT t δ∞∞--∞=-∞⎡⎤=-⎢⎥⎣⎦∑⎰ ()()j e d Ωtan x t t nT t δ∞∞--∞=-∞=-∑⎰()j e ΩnTan x nT ∞-=-∞=∑式中的()a x nT 就是采样后得到的序列()x n , 即()()a x n x nT =()x n 的傅里叶变换为()()j j e enn X x n ωω∞-=-∞=∑比较可知()()j ˆj e aΩTX ΩX ωω==为了在数字计算机上观察分析各种序列的频域特性,通常对()j e X ω在[]0,2π上进行M 点采样来观察分析。
对长度为N 的有限长序列()x n ,有()()1j j 0eekk N nn X x n ωω--==∑其中2π,0,1,,1k k k M Mω==⋅⋅⋅-一个时域离散线性时不变系统的输入/输出关系为()()()()()m y n x n h n x m h n m ∞=-∞=*=-∑上述卷积运算也可以转到频域实现()()()j j j e e e Y X H ωωω= (1.9)三、实验内容及步骤(1) 认真复习采样理论、 离散信号与系统、 线性卷积、 序列的傅里叶变换及性质等有关内容, 阅读本实验原理与方法。
数字信号处理实验报告

数字信号处理实验报告实验一:混叠现象的时域与频域表现实验原理:当采样频率Fs不满足采样定理,会在0.5Fs附近引起频谱混叠,造成频谱分析误差。
实验过程:考虑频率分别为3Hz,7Hz,13Hz 的三个余弦信号,即:g1(t)=cos(6πt), g2(t)=cos(14πt), g3(t)=cos(26πt),当采样频率为10Hz 时,即采样间隔为0.1秒,则产生的序列分别为:g1[n]=cos(0.6πn), g2[n]=cos(1.4πn), g3[n]=cos(2.6πn)对g2[n],g3[n] 稍加变换可得:g2[n]=cos(1.4πn)=cos((2π-0.6π)n)= cos(0.6πn)g3[n]=cos(2.6πn)= cos((2π+0.6π)n)=cos(0.6πn)利用Matlab进行编程:n=1:300;t=(n-1)*1/300;g1=cos(6*pi*t);g2=cos(14*pi*t);g3=cos(26*pi*t);plot(t,g1,t,g2,t,g3);k=1:100;s=k*0.1;q1=cos(6*pi*s);q2=cos(14*pi*s);q3=cos(26*pi*s);hold on; plot(s(1:10),q1(1:10),'bd');figuresubplot(2,2,1);plot(k/10,abs(fft(q1)))subplot(2,2,2);plot(k/10,abs(fft(q2)))subplot(2,2,3);plot(k/10,abs(fft(q3)))通过Matlab软件的图像如图所示:如果将采样频率改为30Hz,则三信号采样后不会发生频率混叠,可运行以下的程序,观察序列的频谱。
程序编程改动如下:k=1:300;q=cos(6*pi*k/30);q1=cos(14*pi*k/30);q2=cos(26*pi*k/30);subplot(2,2,1);plot(k/10,abs(fft(q)))subplot(2,2,2);plot(k/10,abs(fft(q1)))subplot(2,2,3);plot(k/10,abs(fft(q2)))得图像:问题讨论:保证采样后的信号不发生混叠的条件是什么?若信号的最高频率为17Hz,采样频率为30Hz,问是否会发生频率混叠?混叠成频率为多少Hz的信号?编程验证你的想法。
数字信号处理上机实验 作业结果与说明 实验三、四、五

上机频谱分析过程及结果图 上机实验三:IIR 低通数字滤波器的设计姓名:赵晓磊 学号:赵晓磊 班级:02311301 科目:数字信号处理B一、实验目的1、熟悉冲激响应不变法、双线性变换法设计IIR 数字滤波器的方法。
2、观察对实际正弦组合信号的滤波作用。
二、实验内容及要求1、分别编制采用冲激响应不变法、双线性变换法设计巴特沃思、切贝雪夫I 型,切贝雪夫II 型低通IIR 数字滤波器的程序。
要求的指标如下:通带内幅度特性在低于πω3.0=的频率衰减在1dB 内,阻带在πω6.0=到π之间的频率上衰减至少为20dB 。
抽样频率为2KHz ,求出滤波器的单位取样响应,幅频和相频响应,绘出它们的图,并比较滤波性能。
(1)巴特沃斯,双线性变换法Ideal And Designed Lowpass Filter Magnitude Responsefrequency in Hz|H [e x p (j w )]|frequency in pi units|H [ex p (j w )]|Designed Lowpass Filter Phase Response in radians frequency in pi unitsa r g (H [e x p (j w )](2)巴特沃斯,冲激响应不变法(3)切贝雪夫I 型,双线性变换法(4)切贝雪夫Ⅱ型,双线性变换法综合以上实验结果,可以看出,使用不同的模拟滤波器数字化方法时,滤波器的性能可能产生如下差异:使用冲击响应不变法时,使得数字滤波器的冲激响应完全模仿模拟滤波器的冲激响应,也就是时域逼急良好,而且模拟频率和数字频率之间呈线性关系;但频率响应有混叠效应。
frequency in Hz|H [e x p (j w )]|Designed Lowpass Filter Magnitude Response in dBfrequency in pi units|H [e x p (j w )]|frequency in pi unitsa r g (H [e x p (j w )]Ideal And Designed Lowpass Filter Magnitude Responsefrequency in Hz|H [e x p (j w )]|frequency in pi units|H [e xp (j w )]|frequency in pi unitsa r g (H [e x p (j w )]Ideal And Designed Lowpass Filter Magnitude Responsefrequency in Hz|H [e x p (j w )]|frequency in pi units|H [ex p (j w )]|Designed Lowpass Filter Phase Response in radiansfrequency in pi unitsa r g (H [e x p (j w )]使用双线性变换法时,克服了多值映射的关系,避免了频率响应的混叠现象;在零频率附近,频率关系接近于线性关系,高频处有较大的非线性失真。
数字信号处理实验报告五

数字信号处理实验报告实验五:谱分析班级:20110814学号:2011081418姓名:孙明凤日期:2013年12月1日指导教师:田园成绩:实验五谱分析一、实验原理:信号是无限延长的,而在进行信号处理时只能采用有限长信号,所以需要将信号“截断”。
在信号处理中,“截断”被看成是用一个有限长的“窗口”看无限长的信号,或是从分析的角度是无限长的信号乘以有限长的窗函数w(t),有傅里叶变换的性质可知x(t)w(t)↔ 1/2pi *X(jw)*W(jw)如果x(t)是频宽有限信号,而w(t)是频宽无限信号,截断后的信号也必是频宽无锡那信号,从而产生所谓的频谱泄漏。
频谱泄漏是不可避免的但要尽量减少,因此设计了不同的窗函数,满足不同用途的要求,从能量的角度,频谱泄漏也是能量泄漏,因为加窗后,使原来的信号集中在窄频带内的能量分散到无限的频宽范围。
Matlab信号处理工具箱提供了8种窗函数:1)函数boxcar()用于产生矩形窗:w=boxcar(N);2)函数Hanning()用于产生汉宁窗:w=hanning(N);3)函数Hamming()用于产生汉明窗:w=hamming(N);4)函数bartlett()用于产生巴特利窗:w=baetlett(N);5)函数blackman()用于产生布莱克曼窗:w=blackman(N);6)函数tring()用于产生tring窗:w=tring(N);7)函数kaiser()用于产生kaiser窗:w=Kaiser(N);8)函数chebwin用于产生切比雪夫窗:w=chebwin(N);二、实验内容:1、用MA TLAB绘制各种窗函数的形状N=20;w1=boxcar(N);subplot(221);stem(w1);title('boxcar')xlabel('n');ylabel('w1');w2=hanning(N);subplot(222);stem(w2);title('hanning')xlabel('n');ylabel('w2');w3=hamming(N);subplot(223);stem(w3);title('hamming')xlabel('n');ylabel('w3');w4=bartlett(N);subplot(224);stem(w4);title('bartlett')xlabel('n');ylabel('w4');w5=blackman(N);subplot(221);stem(w5);title('blackman')xlabel('n');ylabel('w5');w6=triang(N); subplot(222);stem(w6); title('triang')xlabel('n');ylabel('w6'); beta=60;w7=kaiser(N,beta); subplot(223);stem(w7); title('kaiser')xlabel('n');ylabel('w7'); r=60;w8=chebwin(N,r); subplot(224);stem(w8); title('chebwin')xlabel('n');ylabel('w8');2、用MA TLAB编程绘制各种窗函数的幅频响应N=20;w1=boxcar(N);[H,W]=freqz(w1,1);subplot(221),plot(W/pi/2,abs(H))title('矩形窗振幅特性/dB')xlabel('相对频率');ylabel('|H(W)|');w2=hanning(N);[H,W]=freqz(w2,1);subplot(222),plot(W/pi/2,abs(H))title('汉宁窗振幅特性/dB')xlabel('相对频率');ylabel('|H(W)|');w3=hamming(N);[H,W]=freqz(w3,1);subplot(223),plot(W/pi/2,abs(H))title('汉明窗振幅特性/dB')xlabel('相对频率');ylabel('|H(W)|');w4=bartlett(N);[H,W]=freqz(w4,1);subplot(224),plot(W/pi/2,abs(H))title('巴特利特窗振幅特性/dB')xlabel('相对频率');ylabel('|H(W)|');w5=blackman(N);[H,W]=freqz(w5,1);subplot(221),plot(W/pi/2,abs(H)) title('布莱克曼窗振幅特性/dB') xlabel('相对频率');ylabel('|H(W)|'); w6=triang(N);[H,W]=freqz(w6,1);subplot(222),plot(W/pi/2,abs(H)) title('三角窗振幅特性/dB') xlabel('相对频率');ylabel('|H(W)|'); w7=kaiser(N,60);[H,W]=freqz(w3,1);subplot(223),plot(W/pi/2,abs(H)) title('凯泽窗振幅特性/dB') xlabel('相对频率');ylabel('|H(W)|'); w8=chebwin(N);[H,W]=freqz(w8,1);subplot(224),plot(W/pi/2,abs(H)) title('切比雪夫窗振幅特性/dB') xlabel('相对频率');ylabel('|H(W)|');3、绘制矩形窗的幅频响应,窗长度分别为:N=10,N=20,N=50,N=100。
数字信号处理实验报告(西电)

数字信号处理实验报告班级:****姓名:郭**学号:*****联系方式:*****西安电子科技大学电子工程学院绪论数字信号处理起源于十八世纪的数学,随着信息科学和计算机技术的迅速发展,数字信号处理的理论与应用得到迅速发展,形成一门极其重要的学科。
当今数字信号处理的理论和方法已经得到长足的发展,成为数字化时代的重要支撑,其在各个学科和技术领域中的应用具有悠久的历史,已经渗透到我们生活和工作的各个方面。
数字信号处理相对于模拟信号处理具有许多优点,比如灵活性好,数字信号处理系统的性能取决于系统参数,这些参数很容易修改,并且数字系统可以分时复用,用一套数字系统可以分是处理多路信号;高精度和高稳定性,数字系统的运算字符有足够高的精度,同时数字系统不会随使用环境的变化而变化,尤其使用了超大规模集成的DSP 芯片,简化了设备,更提高了系统稳定性和可靠性;便于开发和升级,由于软件可以方便传送,复制和升级,系统的性能可以得到不断地改善;功能强,数字信号处理不仅能够完成一维信号的处理,还可以试下安多维信号的处理;便于大规模集成,数字部件具有高度的规范性,对电路参数要求不严格,容易大规模集成和生产。
数字信号处理用途广泛,对其进行一系列学习与研究也是非常必要的。
本次通过对几个典型的数字信号实例分析来进一步学习和验证数字信号理论基础。
实验一主要是产生常见的信号序列和对数字信号进行简单处理,如三点滑动平均算法、调幅广播(AM )调制高频正弦信号和线性卷积。
实验二则是通过编程算法来了解DFT 的运算原理以及了解快速傅里叶变换FFT 的方法。
实验三是应用IRR 和FIR 滤波器对实际音频信号进行处理。
实验一●实验目的加深对序列基本知识的掌握理解●实验原理与方法1.几种常见的典型序列:0()1,00,0(){()()(),()sin()j n n n n u n x n Aex n a u n a x n A n σωωϕ+≥<====+单位阶跃序列:复指数序列:实指数序列:为实数 正弦序列:2.序列运算的应用:数字信号处理中经常需要将被加性噪声污染的信号中移除噪声,假定信号 s(n)被噪声d(n)所污染,得到了一个含噪声的信号()()()x n s n d n =+。
数字信号实验报告

数字信号处理上机实验报告姓名:马菲班级:030812 学号:03081161一.实验内容a)利用傅立叶级数展开的方法,自由生成所需的x(t);b)通过选择不同的采样间隔T(分别选T>或<1/2f c),从x(t)获得相应的x(n)(作出x(n)图形);c)对获得的不同x(n)分别作傅立叶变换,分析其频率响应特性(给出幅频与相频特性曲线);d)利用巴特沃思、切比雪夫或椭圆滤波器设计数字滤波器(滤波特性自定),要求通过改变滤波器参数或特性(低通、高通、带通或带阻)设计至少两种数字滤波器,分析所设计滤波器(画出频率特性曲线),并对上述给出的不同x(n)分别进行滤波(画出滤波结果),然后加以讨论;e)利用窗函数设计法或频率采样法设计数字滤波器(滤波特性自定),要求通过改变滤波器参数或特性(低通、高通、带通或带阻等)设计至少两种数字滤波器,分析所设计滤波器(画出频率特性曲线),并对上述给出的不同x(n)分别进行滤波(画出滤波结果),然后加以讨论。
二.实验方法通过MATLAB编程实现生成x(t),以及对其进行上述变换等。
1.生成x(t)利用傅里叶级数展开的方法傅里叶级数:f(t)=A0/2+∑cos(nwn+rn)2通过对N点采样,得到频率响应曲线3利用巴特沃斯设计低通数字滤波器和高通数字滤波器:MATLAB信号处理工具箱函数buttap,buttord和butter是巴特沃斯滤波器设计函数。
4通过汉宁窗和哈明窗设计数字滤波器,调用工具箱函数firl实现窗函数法的设计。
三.部分代码及结果分析1.生成x(t),x(n)及频率响应特性(1)通过傅里叶级数展开方法自定义了一个周期信号y=x(t)如下:%周期信号生成y=3*sin(2*pi*t*10)+0.5*sin(2*pi*t*15)+sin(2*pi*t*40)+sin(2*pi*t*150)(2)选取采样间隔Ts=1400得到相应的序列x(n):%采样序列Ts=1/1400;n=0:279;t=n*Ts;Y=3*sin(2*pi*t*10)+0.5*sin(2*pi*t*15)+sin(2*pi*t*40)+sin(2*pi*t*150); subplot(2,2,2),stem(n,Y);grid on;title('N点采样序列');xlabel('N');ylabel('Y');根据采样间隔取的不同,得到的序列点的密集程度就不一样(3)通过对x(n)作傅里叶变换,得到其幅频特性曲线与相频特性曲线这里以横坐标为采样频率:f=(0:length(Y)-1)'/(length(Y)*Ts);通过傅里叶变换的函数,得到幅频特性曲线:F=abs(fft(Y));相频特性曲线:A=unwrap(angle(fft(Y)));2.通过巴特沃斯设计低通数字滤波器和高通数字滤波器(1)分析:[N,wc]=buttord(wp,ws,rp,rs); 计算阶数N和3dB截止频率。
数字信号处理上机实验报告

数字信号处理上机实验报告实验一熟悉MATLAB环境一、实验目的1、熟悉 MATLAB的主要操作命令。
2、学会简单的矩阵输入和数据读写。
3、掌握简单的绘图命令。
4、用 MATLAB编程并学会创建函数。
5、观察离散系统的频率响应。
二、实验容认真阅读本章附录,在MATLAB环境下重新做一遍附录中的例子,体会各条命令的含义。
在熟悉 MATLAB基本命令的基础上,完成以下实验。
上机实验容:1、数组的加减乘除和乘方运算,输入A1234,B3456,求C A B ,D A B,E A. B,F A./ B,G A.^ B ,并用stem语句画出A、 B、C、 D、 E、F、 G。
程序:>> A=[1 2 3 4];B=[3 4 5 6];C=A+B; D=A-B; E=A.*B; F=A./B; G=A.^B;subplot(2,4,1);stem(A,'.'); subplot(2,4,2);stem(B,'.');subplot(2,4,3);stem(C,'.'); subplot(2,4,4);stem(D,'.');subplot(2,4,5);stem(E,'.'); subplot(2,4,6);stem(F,'.');subplot(2,4,7);stem(G,'.')2、用MATLAB实现下列序列。
a)x(n)0.8n0n15b) x(n)e(0. 2 3 j ) n0n 15c)x(n)3cos(0.125 n0.2 ) 2sin(0.25 n 0.1 ) 0 n 15程序:A)clear;clc;n=[0:15];x1=0.8.^n;subplot(3,1,1),stem(x1)title('x1=0.8^n')xlabel('n'); ylabel('x1');B)clear;clc;n=[0:15];x2=exp((0.2+3j)*n);subplot(3,1,1),stem(x2)title('x2=exp((0.2+3j)*n)')xlabel('n'); ylabel('x2');C)clear;clc;n=[0:15];x3=3*cos(0.125*pi*n+0.2*pi)+2*sin(0.25*pi*n+0.1*pi); subplot(3,1,1),stem(x3)title('x3=3*cos(0.125*pi*n+0.2*pi)+2*sin(0.25*pi*n+0.1*pi)') xlabel('n'); ylabel('x3');3、绘出下列时间常数的图形,对x 轴,y轴以及图形上方均须加上适当的标注:a)x(t )sin( 2 t )0t10sb)x(t )cos(100t )sin(t )0 t 4s>>m=0:0.01:10;n=0:0.01:4;x1t=sin(2*pi*m);x2t=cos(100*pi*n).*sin(pi*n);subplot(2,1,1);plot(m,x1t);subplot(2,1,2);plot(n,x2t);4、给定一因果系统 H(z)=(1+ 2z- 1z-2)/( 1- 0.67z 1z 2),求出并绘制H(z)的幅频响应与相频响应。
西电数字信号处理上机实验报告

数字信号处理上机实验报告14020710021 张吉凯第一次上机实验一:设给定模拟信号()1000t a x t e -=,t 的单位是ms 。
(1) 利用MATLAB 绘制出其时域波形和频谱图(傅里叶变换),估计其等效带宽(忽略谱分量降低到峰值的3%以下的频谱)。
(2) 用两个不同的采样频率对给定的()a x t 进行采样。
○1()()15000s a f x t x n =以样本秒采样得到。
()()11j x n X e ω画出及其频谱。
○2()()11000s a f x t x n =以样本秒采样得到。
()()11j x n X e ω画出及其频谱。
比较两种采样率下的信号频谱,并解释。
(1)MATLAB 程序:N=10; Fs=5; Ts=1/Fs;n=[-N:Ts:N];xn=exp(-abs(n)); w=-4*pi:0.01:4*pi;X=xn*exp(-j*(n'*w));subplot(211)plot(n,xn);title('x_a(t)时域波形');xlabel('t/ms');ylabel('x_a(t)');axis([-10, 10, 0, 1]);subplot(212);plot(w/pi,abs(X));title('x_a(t)频谱图');xlabel('\omega/\pi');ylabel('X_a(e^(j\omega))'); ind = find(X >=0.03*max(X))*0.01;eband = (max(ind) -min(ind));fprintf('等效带宽为%fKHZ\n',eband);运行结果:等效带宽为12.110000KHZ(2)MATLAB程序:N=10;omega=-3*pi:0.01:3*pi;%Fs=5000Fs=5;Ts=1/Fs;n=-N:Ts:N;xn=exp(-abs(n));X=xn*exp(-j*(n'*omega));subplot(2,2,1);stem(n,xn);grid on;axis([-10, 10, 0, 1.25]); title('时域波形(f_s=5000)');xlabel('n');ylabel('x_1(n)');subplot(2,2,2);plot(omega/pi,abs(X));title('频谱图(f_s=5000)');xlabel('\omega/\pi');ylabel('X_1(f)');grid on;%Fs=1000Fs=1;Ts=1/Fs;n=-N:Ts:N;xn=exp(-abs(n));X=xn*exp(-j*(n'*omega));subplot(2,2,3);stem(n,xn);grid on;axis([-10, 10, 0, 1.25]); title('时域波形(f_s=1000)');xlabel('n');ylabel('x_2(n)');grid on; subplot(2,2,4); plot(omega/pi,abs(X)); title('频谱图(f_s=1000)'); xlabel('\omega/\pi'); ylabel('X_2(f)'); grid on;运行结果:实验二:给定一指数型衰减信号()()0cos 2at x t e f t π-=,采样率1s f T=,T 为采样周期。
2017年西电电院数字信号处理上机实验报告六

实验六、FIR数字滤波器设计及其网络结构班级:学号:姓名:成绩:1实验目的(1)熟悉线性相位FIR数字滤波器的时域特点、频域特点和零极点分布;(2)掌握线性相位FIR数字滤波器的窗函数设计法和频率采样设计法;(3)了解IIR数字滤波器和FIR数字滤波器的优缺点及其适用场合。
2 实验内容(1)设计计算机程序.根据滤波器的主要技术指标设计线性相位FIR数字低通、高通、带通和带阻滤波器;(2)绘制滤波器的幅频特性和相频特性曲线.验证滤波器的设计结果是否达到设计指标要求;(3)画出线性相位FIR数字滤波器的网络结构信号流图。
3实验步骤(1)设计相应的四种滤波器的MATLAB程序;(2)画出幅频相频特性曲线;(3)画出信号流图。
4 程序设计%% FIR低通f=[0.2,0.35];m=[1,0];Rp=1;Rs=40;dat1=(10^(Rp/20)-1)/(10^(Rp/20)+1);dat2=10^(-Rs/20);rip=[dat1,dat2];[M,f0,m0,w]=remezord(f,m,rip);M=M+2;hn=remez(M,f0,m0,w);w=0:0.001:pixn=[0:length(hn)-1];H=hn*exp(-j*xn'*w);figuresubplot(2,1,1)plot(w/pi,20*log10(abs(H)));gridon;xlabel('\omega/\pi'),ylabel('|H(e^j^w)|/dB')subplot(2,1,2)plot(w/pi,angle(H)/pi);xlabel('\omega/\pi'),ylabel('\phi(\omega)/\pi')%% FIR高通f=[0.7,0.9];m=[0,1];Rp=1;Rs=60;dat1=(10^(Rp/20)-1)/(10^(Rp/20)+1);dat2=10^(-Rs/20);rip=[dat2,dat1];[M,f0,m0,w]=remezord(f,m,rip);hn=remez(M,f0,m0,w);w=0:0.001:pixn=[0:length(hn)-1];H=hn*exp(-j*xn'*w);figuresubplot(2,1,1)plot(w/pi,20*log10(abs(H)));gridon;xlabel('\omega/\pi'),ylabel('|H(e^j^w)|/dB')subplot(2,1,2)plot(w/pi,angle(H)/pi);xlabel('\omega/\pi'),ylabel('\phi(\omega)/\pi')%% FIR带通f=[0.2,0.35,0.65,0.8];m=[0,1,0];Rp=1;Rs=60;dat1=(10^(Rp/20)-1)/(10^(Rp/20)+1);dat2=10^(-Rs/20);rip=[dat2,dat1,dat2];[M,f0,m0,w]=remezord(f,m,rip);M=M+3hn=remez(M,f0,m0,w);w=0:0.001:pixn=[0:length(hn)-1];H=hn*exp(-j*xn'*w);figuresubplot(2,1,1)plot(w/pi,20*log10(abs(H)));gridon;xlabel('\omega/\pi'),ylabel('|H(e^j^w)|/dB')subplot(2,1,2)plot(w/pi,angle(H)/pi);xlabel('\omega/\pi'),ylabel('\phi(\omega)/\pi')%% FIR带阻f=[0.2,0.35,0.65,0.8];m=[1,0,1];Rp=1;Rs=60;dat1=(10^(Rp/20)-1)/(10^(Rp/20)+1);dat2=10^(-Rs/20);rip=[dat1,dat2,dat1];[M,f0,m0,w]=remezord(f,m,rip);hn=remez(M,f0,m0,w);w=0:0.001:pixn=[0:length(hn)-1];H=hn*exp(-j*xn'*w);figuresubplot(2,1,1)plot(w/pi,20*log10(abs(H)));gridon;xlabel('\omega/\pi'),ylabel('|H(e^j^w)|/dB')subplot(2,1,2)plot(w/pi,angle(H)/pi);xlabel('\omega/\pi'),ylabel('\phi(\omega)/\pi') 5实验结果及分析(1)FIR低通滤波器自动得到的M值不满足要求.故我们将M加上2 在w=0.2π时.H=-0.5dB;w=0.35π时.H=-41dB。
西电数字信号处理上机实验报告

数字信号处理上机实验报告14020710021 张吉凯第一次上机实验一:设给定模拟信号()1000t a x t e -=,t 的单位是ms 。
(1) 利用MATLAB 绘制出其时域波形和频谱图(傅里叶变换),估计其等效带宽(忽略谱分量降低到峰值的3%以下的频谱)。
(2) 用两个不同的采样频率对给定的()a x t 进行采样。
○1()()15000s a f x t x n =以样本秒采样得到。
()()11j x n X e ω画出及其频谱。
○2()()11000s a f x t x n =以样本秒采样得到。
()()11j x n X e ω画出及其频谱。
比较两种采样率下的信号频谱,并解释。
(1)MATLAB 程序:N=10; Fs=5; T s=1/Fs;n=[-N:T s:N];xn=exp(-abs(n)); w=-4*pi:0.01:4*pi; X=xn*exp(-j*(n'*w)); subplot(211) plot(n,xn);title('x_a(t)时域波形');xlabel('t/ms');ylabel('x_a(t)'); axis([-10, 10, 0, 1]); subplot(212);plot(w/pi,abs(X)); title('x_a(t)频谱图');xlabel('\omega/\pi');ylabel('X_a(e^(j\omega))');ind = find(X >=0.03*max(X))*0.01; eband = (max(ind) -min(ind)); fprintf('等效带宽为%fKHZ\n',eband); 运行结果:等效带宽为12.110000KHZ(2)MATLAB程序:N=10;omega=-3*pi:0.01:3*pi;%Fs=5000Fs=5;T s=1/Fs;n=-N:T s:N;xn=exp(-abs(n));X=xn*exp(-j*(n'*omega));subplot(2,2,1);stem(n,xn);grid on;axis([-10, 10, 0, 1.25]); title('时域波形(f_s=5000)');xlabel('n');ylabel('x_1(n)');subplot(2,2,2);plot(omega/pi,abs(X));title('频谱图(f_s=5000)');xlabel('\omega/\pi');ylabel('X_1(f)');grid on;%Fs=1000Fs=1;T s=1/Fs;n=-N:T s:N;xn=exp(-abs(n));X=xn*exp(-j*(n'*omega));subplot(2,2,3);stem(n,xn);grid on;axis([-10, 10, 0, 1.25]); title('时域波形(f_s=1000)');xlabel('n');ylabel('x_2(n)');grid on;subplot(2,2,4);plot(omega/pi,abs(X));title('频谱图(f_s=1000)');xlabel('\omega/\pi');ylabel('X_2(f)');grid on;运行结果:实验二:给定一指数型衰减信号()()0cos 2at x t e f t π-=,采样率1s f T=,T 为采样周期。
西电数字信号处理上机实验

实验一1-1、a=[-2 0 1 -1 3];b=[1 2 0 -1];c=conv(a,b);M=length(c)-1;n=0:1:M;stem(n,c);xlabel('n');ylabel('幅度');title('离散卷积’);1-2、N=41;a=[0.8 -0.44 0.36 0.22]; b=[1 0.7 -0.45 -0.6];x=[1 zeros(1,N-1)];k=0:1:N-1;y=filter(a,b,x);stem(k,y)xlabel('n');ylabel('幅度'); title('差分方程');1-3、k=256;num=[0.8 -0.44 0.36 0.02];den=[1 0.7 -0.45 -0.6];w=0:pi/k:pi;h=freqz(num,den,w);subplot(2,2,1);plot(w/pi,real(h));gridtitle('实部');xlabel('\omega/\pi');ylabel('幅度'); subplot(2,2,2);plot(w/pi,imag(h));gridtitle('虚部');xlabel('\omega/\pi');ylabel('Amplitude'); subplot(2,2,3);plot(w/pi,abs(h));gridtitle('幅度谱');xlabel('\omega/\pi');ylabel('幅值'); subplot(2,2,4);plot(w/pi,angle(h));gridtitle('相位谱');xlabel('\omega/\pi');ylabel('弧度');实验二2-1、N=16;n=0:1:15;p=8;q=4;a=0.1;f=0.0625;xa=exp(-((n-p).^2)./q);figure(1)stem(n, xa,'.');title('xa(n)序列')xlabel('n')ylabel('xa(n)')grid on[H, w] = freqz(xa, 1, [], 'whole', 1); Hamplitude = abs(H);Hphase = angle(H);Hphase = unwrap(Hphase);figure(2)subplot(2, 1, 1)plot(w, Hamplitude)title('幅频响应')xlabel('w/(2*pi)')ylabel('|H(exp(jw))|') grid onsubplot(2, 1, 2)plot(w, Hphase)title('相频响应')xlabel('w/(2*pi)')ylabel('fai(H(exp(jw)))') grid on2-2、n=0:1:15;a=0.1;f1=0.0625;f2=0.04375;f3=0.05625;xb1=exp(-a*n).*sin(2*pi*f1*n);figuresubplot(3,2,1)stem(n, xb1,'.');title('f=0.0625的时域特性')xlabel('n')ylabel('xb1(n)')grid on[H, w] = freqz(xb1, 1, [], 'whole', 1); Hamplitude = abs(H);subplot(3,2,2)plot(w, Hamplitude)title('f=0.0625的幅频响应')xlabel('w/(2*pi)')ylabel('|H(exp(jw))|')grid onxb2=exp(-a*n).*sin(2*pi*f2*n);subplot(3,2,3)stem(n, xb2,'.');title('f=0.04375的时域特性')xlabel('n')ylabel('xb2(n)')grid on[H, w] = freqz(xb2, 1, [], 'whole', 1); Hamplitude = abs(H);subplot(3,2,4)plot(w, Hamplitude)title('f=0.04375的幅频响应')xlabel('w/(2*pi)')ylabel('|H(exp(jw))|')grid onxb3=exp(-a*n).*sin(2*pi*f3*n);subplot(3,2,5)stem(n, xb3,'.');title('f=0.05625的时域特性')xlabel('n')ylabel('xb3(n)')grid on[H, w] = freqz(xb3, 1, [], 'whole', 1); Hamplitude = abs(H);subplot(3,2,6)plot(w, Hamplitude)title('f=0.05625的幅频响应')xlabel('w/(2*pi)')ylabel('|H(exp(jw))|')grid on2-3、n1=0:1:3;xc1=n1+1;n2=4:7;xc2=8-n2;xc=[xc1,xc2];n =[n1,n2];figurestem(n,xc);xlabel('n'); ylabel('xc');title('三角序列');n1=0:1:3;xd1=4-n1;n2=4:7;xd2=n2-3;xd=[xd1,xd2];n =[n1,n2];figurestem(n,xd);xlabel('n'); ylabel('xd');title('反三角序列');N = 16;[H1,w1] = freqz(xc,1, 256, 'whole', 1); Hamplitude1 = abs(H1);figureplot(2*w1, Hamplitude1)title('xc幅频响应')xlabel('w/pi')ylabel('|H(exp(jw))|')grid on[H2,w2] = freqz(xd,1, 256, 'whole', 1); Hamplitude2 = abs(H2);figureplot(2*w2, Hamplitude2)title('xd幅频响应')xlabel('w/pi')ylabel('|H(exp(jw))|')grid on[H3, w3] = freqz(xc, 1, N, 'whole', 1); Hamplitude3 = abs(H3);figuresubplot(2, 1, 1)h3 = stem(2*w3, Hamplitude3, '*');title('xc幅频响应进行N点FFT’);xlabel('n')ylabel('|H(exp(jw))|')grid on[H4, w4] = freqz(xd, 1, N, 'whole', 1); Hamplitude4 = abs(H4);subplot(2, 1, 2)h4 = stem(2*w4, Hamplitude4, '*');title('xd幅频响应进行N点FFT');xlabel('n')ylabel('|H(exp(jw))|')grid on2-4、N = 128;f1 = 1/16;n = 0:N-1;xn = sin(2*pi*0.125.*n)+ cos(2*pi*(0.125+f1).*n); figurestem(n,xn);figuresubplot(2,1,1),plot(n,abs(fft(xn)));title('f =1/16 幅频响应');f2 = 1/64;xn = sin(2*pi*0.125.*n)+ cos(2*pi*(0.125+f2).*n); subplot(2,1,2),plot(n,abs(fft(xn)));title('f =1/64 幅频响应');2-5、N=16;n=0:1:15;p=8;q=2;a=0.1;f=0.0625;xa=exp(-((n-p).^2)./q);xb=exp(-a*n).*sin(2*pi*f*n);%线性卷积x=conv(xa,xb);XDft= fft(x, 32);XDftR = abs(XDft);XDftPhase = angle(XDft);XDftPhase = unwrap(XDftPhase);figure(1);stem(x,'.');title('x(n)序列');xlabel('n')ylabel('x(n)')grid onfigure(2)subplot(2, 1, 1)stem(XDftR, '.');title('X(k)的幅度’);xlabel('k')ylabel('|X(k)|')grid onsubplot(2, 1, 2)stem(XDftPhase, '.');title('X(k)的相角')xlabel('k')ylabel('fai((X(k)))')grid on%圆周卷积XDft161 = fft(xa, N);XDft16R1 = abs(XDft161);XDft16Phase1 = angle(XDft161);XDft16Phase1 = unwrap(XDft16Phase1); XDft162 = fft(xb, N);XDft16R2 = abs(XDft162);XDft16Phase2 = angle(XDft162);XDft16Phase2 = unwrap(XDft16Phase2); XDft16=XDft161.*XDft162;XDft16R=XDft16R1.*XDft16R2;XDft16Phase=XDft16Phase2 +XDft16Phase1 ; x = ifft(XDft16, N);figure(3)stem(x,'.')title('x(n)序列')xlabel('n')ylabel('x(n)')grid onfigure(4)subplot(2, 1, 1)t= 0 : 1 : N - 1;stem(t, XDft16R, '.');title('X(k)的幅度')xlabel('k')ylabel('|X(k)|')grid onsubplot(2, 1, 2)stem(t,XDft16Phase, '.');title('X(k)的相角')xlabel('k')ylabel('fai((X(k)))')grid on2-6、xe=rand(1,512);n1=0:1:3;xc1=n1+1;n2=4:7;xc2=8-n2;xc=[xc1,xc2];%重叠相加法yn=zeros(1,519);for j=0:7xj=xe(64*j+1:64*(j+1));xak=fft(xj,71);xck=fft(xc,71);yn1=ifft(xak.*xck);temp=zeros(1,519);temp(64*j+1:64*j+71)=yn1; yn=yn+temp;end;n=0:518;figure(1)subplot(2,1,1);plot(n,yn);xlabel('n');ylabel('y(n)');title('xc(n)与xe(n)的线性卷积的时域波形-重叠相加法'); subplot(2,1,2);plot(n,abs(fft(yn)));xlabel('k');ylabel('Y(k)');axis([0,600,0,300]);title('xc(n)Óëxe(n)的线性卷积的幅频特性-重叠相加法'); %重叠保留法k=1:7;xe1=k-k;xe_1=[xe1,xe];yn_1=zeros(1,519);for j=0:7xj_1=xe_1(64*j+1:64*j+71);xak_1=fft(xj_1);xck_1=fft(xc,71);yn1_1=ifft(xak_1.*xck_1);temp_1=zeros(1,519);temp_1(64*j+1:64*j+64)=yn1_1(8:71);yn_1=yn_1+temp_1;end;n=0:518;figure(2)subplot(2,1,1);plot(n,yn_1);xlabel('n');ylabel('y(n)');title(' xc(n)的线性卷积的时域波形-重叠保留法'); subplot(2,1,2);plot(n,abs(fft(yn_1)));xlabel('k');ylabel('Y(k)');axis([0,600,0,300]);title('xc(n)Óëxe(n)的线性卷积的幅频特性-重叠保留法');实验三3-1、Wp=0.3;Ws=0.2;Rp=0.8;Rs=20;[N,Wpo]=cheb1ord(Wp,Ws,Rp,Rs);[Bz,Az]=cheby1(N,Rp,Wpo,'high');w=0:0.1:pi;[H,w1]=freqz(Bz,Az,w);H=20*log10(abs(H));plot(w/pi,H),grid onxlabel('\omega/\pi');ylabel('|H(e^j^\omega)|/dB') title('Chebyshev高通滤波器');3-2、Wp=0.2;Ws=0.3;Rp=1;Rs=25;[N,Wc]=buttord(Wp,Ws,Rp,Rs);[Bs,As]=butter(N,Wc,'s');[Bz,Az]=impinvar(Bs,As);w=0:0.1:pi;[H,w1]=freqz(Bz,Az,w);H=20*log10(abs(H));subplot(211);plot(w/pi,H),grid onxlabel('\omega/\pi');ylabel('|H(e^j^\omega)|/dB') title('脉冲响应不变法')[N,Wc]=buttord(Wp,Ws,Rp,Rs);[Bz,Az]=butter(N,Wc);w=0:0.1:pi;[H,w1]=freqz(Bz,Az,w);H=20*log10(abs(H));subplot(212);plot(w/pi,H),grid onxlabel('\omega/\pi');ylabel('|H(e^j^\omega)|/dB') title('双线性变换法')3-3、Wp=1.2/8;Ws=2/8;Rp=0.5;Rs=40;[N,Wpo]=cheb1ord(Wp,Ws,Rp,Rs);[Bz,Az]=cheby1(N,Rp,Wpo);w=0:0.1:pi;[H,w1]=freqz(Bz,Az,w);H=20*log10(abs(H));subplot(311);plot(w/pi,H),grid onxlabel('\omega/\pi');ylabel('|H(e^j^\omega)|/dB') title('切比雪夫')[N,Wc]=buttord(Wp,Ws,Rp,Rs);[Bz,Az]=butter(N,Wc);w=0:0.1:pi;[H,w1]=freqz(Bz,Az,w);H=20*log10(abs(H));subplot(312);plot(w/pi,H),grid onxlabel('\omega/\pi');ylabel('|H(e^j^\omega)|/dB') title('巴特沃斯')[N,Wpo]=ellipord(Wp,Ws,Rp,Rs);[Bz,Az]=ellip(N,Rp,Rs,Wpo);w=0:0.1:pi;[H,w1]=freqz(Bz,Az,w);H=20*log10(abs(H));subplot(313);plot(w/pi,H),grid ontitle('椭圆')xlabel('\omega/\pi');ylabel('|H(e^j^\omega)|/dB')3-4、Wp1=2/15;Wpu=0.2;Ws1=0.1;Wsu=0.4;Rp=3;Rs=20;Wp=[Wp1,Wpu];Ws=[Ws1,Wsu];[N,Wc]=buttord(Wp,Ws,Rp,Rs);[Bz,Az]=butter(N,Wc);w=0:0.1:pi;[H,w1]=freqz(Bz,Az,w);H=20*log10(abs(H));plot(w/pi,H),grid onxlabel('\omega/\pi');ylabel('|H(e^j^\omega)|/dB') title('双线性变换法Butterworth型数字带通滤波器')。
西电数字信号处理大作业

实验一、信号的采样clc,clear;dt=0.001;tf=6;t=0:dt:tf;xa=sqrt(t)+cos(t);T=0.5;n=0:tf/T;x=sqrt(n*T)+cos(n*T);figure(1)subplot(2,1,1)plot(t,xa),grid on ;title('original image')subplot(2,1,2)stem(n*T,x),grid on ,title('digital image')实验二、信号与系统的时域分析差分方程为)()2()1()(21n bx n y a n y a n y +----=,其中8.01-=a ,64.02=a ,866.0=b 。
系统单位脉冲响应)(n ha1=-0.8;a2=0.64;b=0.866;ys=0;xn=[1,zeros(1,49)];B=1;A=[1,a1,a2];xi=filtic(B,A,ys);yn=filter(B,A,xn,xi);n=0:length(yn)-1;subplot(1,1,1);stem(n,yn,'.')title('(a)');xlabel('n');ylabel('y(n)')输入x(n)=cos(n)T=0.1;z=cos(n*T);zn=conv(yn,z); figure(2);n1=1:99;stem(n1,zn,'.')实验三、系统的频域和Z域分析程序代码(画出dtft的幅度和频率谱)clc,clear;n=0:1:7;x=(0.9*exp(j*pi/3)).^n;w=0:pi/200:pi;X=x*exp(-j).^(n'*w);realX=real(X);imagX=imag(X);angX=angle(X);magX=abs(X);subplot(2,2,1);plot(w/pi,magX);grid xlabel('frequency in pi unit');title('magnitude part');subplot(2,2,2);plot(w/pi,realX);grid xlabel('frequency in pi unit');title('real part');subplot(2,2,3);plot(w/pi,imagX);grid xlabel('frequency in pi unit');title('imaginary part');subplot(2,2,4);plot(w/pi,angX);grid xlabel('frequency in pi unit');title('angel part');clc,clear;a=[1,-0.5,0.06];b=[1,1,0];m=0:length(b)-1;l=0:length(a)-1;w=0:pi/500:pi;num=b*exp(-j*m'*w);den=a*exp(-j*l'*w);H=num./den;magH=abs(H);angH=angle(H);H1=freqz(b,a,w);magH1=abs(H1);angH1=angle(H1);subplot(2,2,2);plot(w/pi,angH);grid;xlabel('w(frequency in pi units)');ylabel('Ïàλrad/w');subplot(2,2,1);plot(w/pi,magH);grid;xlabel('w(frequency in pi units)');ylabel('·ù¶È|H|');subplot(2,2,3);plot(w/pi,magH1);grid;xlabel('w(frequency in pi units)');ylabel('·ù¶È|H1|');subplot(2,2,4);plot(w/pi,angH);grid;xlabel('w(frequency in pi units)');ylabel('Ïàλrad/w');axis([0,1,-0.8,0]); figure(2);zplane(b,a);实验四、信号的频谱分析程序代码clc,clear;n=0:7;k=0:7;N=8;w=n*(2*pi)/8;x=(0.9*exp(j*pi/3)).^n;X1=[x zeros(1,8)];X2=[X1 zeros(1,16)];XK=x*exp(-j*k'*w);k1=0:15;n1=0:15;w1=n1*(2*pi)/16;XK1=X1*exp(-j*k1'*w1);k2=0:31;n2=0:31;w2=n2*(2*pi)/16;XK2=X2*exp(-j*k2'*w2);w3=0:pi/200:2*pi;X=x*exp(-j*n'*w3);magX=abs(X);angX=angle(X);magXK=abs(XK);angXK=angle(XK);magXK1=abs(XK1);angXK1=angle(XK1);magXK2=abs(XK2);angXK2=angle(XK2);subplot(4,2,1);plot(w3/pi,magX);xlabel('w/pi');ylabel('·ù¶È|X|');grid on;subplot(4,2,2);plot(w3/pi,angX);xlabel('w/pi');ylabel('Ïàλrad/pi'); subplot(4,2,3);stem(n,magXK);xlabel('K');ylabel('·ù¶È|XK|');subplot(4,2,4);stem(n,magXK);xlabel('K');ylabel('Ïàλrad/pi'); subplot(4,2,5);stem(n1,magXK1);xlabel('K1');ylabel('·ù¶È|XK1|'); subplot(4,2,6);stem(n1,magXK1);xlabel('K1');ylabel('Ïàλrad/pi'); subplot(4,2,7);stem(n2,magXK2);xlabel('K2');ylabel('·ù¶È|XK2|'); subplot(4,2,8);stem(n2,magXK2);xlabel('K2');ylabel('Ïàλrad/pi');实验五、IIR数字滤波器设计IIR汉宁窗低通高通低通巴特沃斯通带截止频率wp=0.2pi 通带最大衰减R=1dB阻带截止频率wp=0.35pi 阻带最大衰减R=10dBclc,clear;Wp=0.2;Ws=0.35;Rp=1;Rs=100;[N,Wc]=buttord(Wp,Ws,Rp,Rs);[Bz,Az]=butter(N,Wc)w=0:0.1:pi;[H,w1]=freqz(Bz,Az,w);;ang=angle(H);H=20*log10(abs(H))subplot(4,2,1); plot(w/pi,H) ;gridon ;xlabel('\omega/\pi');ylabel('|H(e^j^\omega)|/dB');title('µÍͨÂ˲¨Æ÷')subplot(4,2,2);plot(w/pi,ang);gridon ;xlabel('\omega/\pi');ylabel('Phase/dB')[Bz1,Az1]=butter(N,Wc,'high')w=0:0.1:pi;[H1,w2]=freqz(Bz1,Az1,w);ang1=angle(H1);H1=20*log10(abs(H1))subplot(4,2,3); plot(w/pi,H1) ;gridon ;xlabel('\omega/\pi');ylabel('|H(e^j^\omega)|/dB');title('¸ßͨÂ˲¨Æ÷')subplot(4,2,4);plot(w/pi,ang1);gridon ;xlabel('\omega/\pi');ylabel('Phase/dB')Wp1=[0.2 0.8];Ws1=[0.35 0.65];[N2,Wc1]=buttord(Wp1,Ws1,Rp,Rs);[Bz2,Az2]=butter(N2,Wc1,'stop')w=0:0.1:pi;[H2,w3]=freqz(Bz2,Az2,w);ang2=angle(H2);H2=20*log10(abs(H2))subplot(4,2,5); plot(w/pi,H2) ;gridon ;xlabel('\omega/\pi');ylabel('|H(e^j^\omega)|/dB');title('´ø×èÂ˲¨Æ÷')subplot(4,2,6);plot(w/pi,ang2);gridon ;xlabel('\omega/\pi');ylabel('Phase/dB')Wp1=[0.2 0.8];Ws1=[0.35 0.65];[N2,Wc1]=buttord(Wp1,Ws1,Rp,Rs);[Bz3,Az3]=butter(N2,Wc1)w=0:0.1:pi;[H3,w4]=freqz(Bz3,Az3,w);ang3=angle(H3);H3=20*log10(abs(H3))subplot(4,2,7); plot(w/pi,H3) ;gridon ;xlabel('\omega/\pi');ylabel('|H(e^j^\omega)|/dB');title('´øͨÂ˲¨Æ÷')subplot(4,2,8);plot(w/pi,ang3);gridon ;xlabel('\omega/\pi');ylabel('Phase/dB')切比雪夫1型通带截止频率wp=0.7pi 通带最大衰减R=1dB阻带截止频率wp=0.5pi 阻带最大衰减R=40dBclc,clear;Wp=0.7;Ws=0.5;Rp=1;Rs=40;[N,Wpo]=cheb1ord(Wp,Ws,Rp,Rs);[Bz,Az]=cheby1(N,Rp,Wpo)w=0:0.1:pi;[H,w1]=freqz(Bz,Az,w);ang=angle(H);H=20*log10(abs(H))subplot(4,2,1); plot(w/pi,H) ;gridon ;xlabel('\omega/\pi');ylabel('|H(e^j^\omega)|/dB');title('µÍͨÂ˲¨Æ÷')subplot(4,2,2);plot(w/pi,ang);gridon ;xlabel('\omega/\pi');ylabel('Phase/dB')[Bz1,Az1]=cheby1(N,Rp,Wpo,'high');w=0:0.1:pi;[H1,w2]=freqz(Bz1,Az1,w);ang1=angle(H1);H1=20*log10(abs(H1))subplot(4,2,3); plot(w/pi,H1) ;gridon ;xlabel('\omega/\pi');ylabel('|H(e^j^\omega)|/dB');title('¸ßͨÂ˲¨Æ÷')subplot(4,2,4);plot(w/pi,ang1);gridon ;xlabel('\omega/\pi');ylabel('Phase/dB')Wp1=[0.2 0.8];Ws1=[0.35 0.65];[N2,Wpo1]=cheb1ord(Wp1,Ws1,Rp,Rs);[Bz2,Az2]=cheby1(N2,Rp,Wpo1,'stop')w=0:0.1:pi;[H2,w3]=freqz(Bz2,Az2,w);ang2=angle(H2);H2=20*log10(abs(H2))subplot(4,2,5); plot(w/pi,H2) ;gridon ;xlabel('\omega/\pi');ylabel('|H(e^j^\omega)|/dB');title('´ø×èÂ˲¨Æ÷')subplot(4,2,6);plot(w/pi,ang2);gridon ;xlabel('\omega/\pi');ylabel('Phase/dB')Wp1=[0.2 0.8];Ws1=[0.35 0.65];[N2,Wpo1]=cheb1ord(Wp1,Ws1,Rp,Rs);[Bz3,Az3]=cheby1(N2,Rp,Wpo1)w=0:0.1:pi;[H3,w4]=freqz(Bz3,Az3,w);ang3=angle(H3);H3=20*log10(abs(H3))subplot(4,2,7); plot(w/pi,H3) ;gridon ;xlabel('\omega/\pi');ylabel('|H(e^j^\omega)|/dB');title('´øͨÂ˲¨Æ÷')subplot(4,2,8);plot(w/pi,ang3);gridon ;xlabel('\omega/\pi');ylabel('Phase/dB')实验六、FIR数字滤波器设计FIR汉宁窗低通高通低通% 采用Hamming窗设计一个带阻FIR滤波器阻带:0~0.5pi,阻带最小衰减Rs=40dB;通带:0.5~pi,通带最大衰减:Rp=1dB。
西安电子科技大学数字信号处理上机报告

数字信号处理上机实验报告班级:020915**:***学号:********实验一:信号、系统及系统响应1、实验目的(1) 熟悉连续信号经理想采样前后的频谱变化关系,加深对时域采样定理的理解。
(2) 熟悉时域离散系统的时域特性。
(3) 利用卷积方法观察分析系统的时域特性。
(4) 掌握序列傅里叶变换的计算机实现方法,利用序列的傅里叶变换对连续信号、离散信号及系统响应进行频域分析。
2、实验原理与方法(1) 时域采样。
(2) LTI系统的输入输出关系。
3、实验内容及步骤(1) 认真复习采样理论、离散信号与系统、线性卷积、序列的傅里叶变换及性质等有关内容,阅读本实验原理与方法。
(2) 编制实验用主程序及相应子程序。
①信号产生子程序,用于产生实验中要用到的下列信号序列:a. xa(t)=A*e^-at *sin(Ω0t)u(t)b. 单位脉冲序列:xb(n)=δ(n)c. 矩形序列:xc(n)=RN(n), N=10②系统单位脉冲响应序列产生子程序。
本实验要用到两种FIR系统。
a. h a(n)=R10(n);b. h b(n)=δ(n)+2.5δ(n-1)+2.5δ(n-2)+δ(n-3)③有限长序列线性卷积子程序用于完成两个给定长度的序列的卷积。
可以直接调用MATLAB语言中的卷积函数conv。
conv用于两个有限长度序列的卷积,它假定两个序列都从n=0 开始。
调用格式如下:y=conv (x, h)4、实验结果记录①分析采样序列的特性。
a. 取采样频率fs=1 kHz,,即T=1 ms。
b. 改变采样频率,fs=300 Hz,观察|X(e^jω)|的变化,并做记录(打印曲线);进一步降低采样频率,fs=200 Hz,观察频谱混叠是否明显存在,说明原因,并记录(打印)这时的|X(e^jω)|曲线。
程序代码如下:A=444.128;a=50*sqrt(2)*pi;m=50*sqrt(2)*pi;fs1=1000; %此时,取采样频率fs=1kHz,之后,改变采样频率分别为200Hz,300Hz fs2=300;fs3=200;T1=1/fs1;T2=1/fs2;T3=1/fs3;N=100;n=[0:N-1];x1=A*exp(-a*n*T1).*sin(m*n*T1);x2=A*exp(-a*n*T2).*sin(m*n*T2);x3=A*exp(-a*n*T3).*sin(m*n*T3);w=linspace(-pi,pi,10000); %设置w的范围X1=x1*exp(-j*n'*w); %对x1(n)做DTFT变换X2=x2*exp(-j*n'*w); %对x2(n)做DTFT变换X3=x3*exp(-j*n'*w); %对x3(n)做DTFT变换figure(1)subplot(1,3,1)plot(w/pi,abs(X1)); %绘制x1(n)的幅度谱xlabel('\omega/π');ylabel('|H(e^j^\omega)|')title('采样频率为1000Hz时的幅度谱');subplot(1,3,2)plot(w/pi,abs(X2)); %绘制x2(n)的幅度谱xlabel('\omega/π');ylabel('|H(e^j^\omega)|')title('采样频率为300Hz时的幅度谱');subplot(1,3,3)plot(w/pi,abs(X3)); %绘制x3(n)的幅度谱xlabel('\omega/π');ylabel('|H(e^j^\omega)|')title('采样频率为200Hz时的幅度谱');②时域离散信号、系统和系统响应分析。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
实验五、IIR数字滤波器设计及其网络结构
班级:学号:姓名:成绩:
1实验目的
(1)熟悉数字滤波的基本概念、数字滤波器的主要技术指标及其物理意义;
(2)掌握巴特沃斯和切比雪夫模拟低通滤波器的设计方法和IIR数字低通滤波器的脉冲响应不变设计法、双线性变换法设计方法。
(3)了解模拟和数字滤波器的频率变换、IIR数字滤波器的直接(优化)设计方法;
2 实验内容
(1)设计计算机程序,根据滤波器的主要技术指标设计IIR数字巴特沃斯和切比雪夫低通、高通、带通和带阻滤波器;
(2)绘制滤波器的幅频特性和相频特性曲线,验证滤波器的设计结果是否达到设计指标要求;
(3)画出数字滤波器的直接型、级联型、并联型网络结构信号流图。
3实验步骤
(1)设计相应的八种滤波器的MATLAB程序;
(2)画出幅频相频特性曲线;
(3)画出信号流图。
4程序设计
%% 巴特沃斯低通
wp=0.2;ws=0.35;rp=1;rs=10;
[N,wc]=buttord(wp,ws,rp,rs);
[B,A]=butter(N,wc);
w=0:0.001:pi;
[H,w]=freqz(B,A,w);
H1=20*log10(abs(H))
subplot(2,1,1)
plot(w/pi,H1),grid on;xlabel('\omega/\pi'),ylabel('|H(e^i^\omega)|')
subplot(2,1,2)
plot(w/pi,angle(H)/pi);xlabel('\omega/\pi'),ylabel('\phi(\omega)/\pi')
%% 巴特沃斯高通
wp=0.8;ws=0.6;rp=1;rs=10;
[N,wc]=buttord(wp,ws,rp,rs);
[B,A]=butter(N,wc,'high');
w=0:0.001:pi;
[H,w]=freqz(B,A,w);
H1=20*log10(abs(H));
subplot(2,1,1)
plot(w/pi,H1),grid on;
xlabel('\omega/\pi'),ylabel('|H(e^i^\omega)|')
subplot(2,1,2)
plot(w/pi,angle(H)/pi);xlabel('\omega/\pi'),ylabel('\phi(\omega)/\pi') %% 巴特沃斯带通
wpl=0.4;wpu=0.6;wsl=0.2;wsu=0.8
wp=[wpl,wpu];ws=[wsl,wsu];rp=1;rs=20;
[N,wc]=buttord(wp,ws,rp,rs);
[B,A]=butter(N,wc);
w=0:0.001:pi;
[H,w]=freqz(B,A,w);
H1=20*log10(abs(H));
subplot(2,1,1)
plot(w/pi,H1),grid on;xlabel('\omega/\pi'),ylabel('|H(e^i^\omega)|') subplot(2,1,2)
plot(w/pi,angle(H)/pi);xlabel('\omega/\pi'),ylabel('\phi(\omega)/\pi') %% 巴特沃斯带阻
wpl=0.2;wpu=0.8;wsl=0.4;wsu=0.6
wp=[wpl,wpu];ws=[wsl,wsu];rp=1;rs=20;
[N,wc]=buttord(wp,ws,rp,rs);
[B,A]=butter(N,wc,'stop');
w=0:0.001:pi;
[H,w]=freqz(B,A,w);
H1=20*log10(abs(H));
subplot(2,1,1)
plot(w/pi,H1),grid on;
xlabel('\omega/\pi'),ylabel('|H(e^i^\omega)|')
subplot(2,1,2)
plot(w/pi,angle(H)/pi);xlabel('\omega/\pi'),ylabel('\phi(\omega)/\pi') %% 切比雪夫低通
wp=0.2;ws=0.5;rp=1;rs=40;
[N,wpo]=cheb1ord(wp,ws,rp,rs);
[B,A]=cheby1(N,rp,wpo);
w=0:0.001:pi;
[H,w]=freqz(B,A,w);
H1=20*log10(abs(H));
subplot(2,1,1)
plot(w/pi,H1),grid on;xlabel('\omega/\pi'),ylabel('|H(e^i^\omega)|')
subplot(2,1,2)
plot(w/pi,angle(H)/pi);xlabel('\omega/\pi'),ylabel('\phi(\omega)/\pi')
%% 切比雪夫高通
wp=0.7;ws=0.5;rp=1;rs=40;
[N,wpo]=cheb1ord(wp,ws,rp,rs);
[B,A]=cheby1(N,rp,wpo,'high');
w=0:0.001:pi;
[H,w]=freqz(B,A,w);
H1=20*log10(abs(H));
subplot(2,1,1)
plot(w/pi,H1),grid on;xlabel('\omega/\pi'),ylabel('|H(e^i^\omega)|')
subplot(2,1,2)
plot(w/pi,angle(H)/pi);xlabel('\omega/\pi'),ylabel('\phi(\omega)/\pi')
%% 切比雪夫带通
wpl=0.4;wpu=0.6;wsl=0.2;wsu=0.8
wp=[wpl,wpu];ws=[wsl,wsu];rp=1;rs=20;
[N,wpo]=cheb1ord(wp,ws,rp,rs);
[B,A]=cheby1(N,rp,wpo);
w=0:0.001:pi;
[H,w]=freqz(B,A,w);
H1=20*log10(abs(H));
subplot(2,1,1)
plot(w/pi,H1),grid on;xlabel('\omega/\pi'),ylabel('|H(e^i^\omega)|')
subplot(2,1,2)
plot(w/pi,angle(H)/pi);xlabel('\omega/\pi'),ylabel('\phi(\omega)/\pi')
%% 切比雪夫带阻
wpl=0.2;wpu=0.8;wsl=0.4;wsu=0.6
wp=[wpl,wpu];ws=[wsl,wsu];rp=1;rs=20;
[N,wpo]=cheb1ord(wp,ws,rp,rs);
[B,A]=cheby1(N,rp,wpo,'stop');
w=0:0.001:pi;
[H,w]=freqz(B,A,w);
H1=20*log10(abs(H));
subplot(2,1,1)
plot(w/pi,H1),grid on;xlabel('\omega/\pi'),ylabel('|H(e^i^\omega)|')
subplot(2,1,2)
plot(w/pi,angle(H)/pi);xlabel('\omega/\pi'),ylabel('\phi(\omega)/\pi')
5实验结果及分析
(1)巴特沃斯低通
W=0.5πi时,H=-0.75dB,w=0.35π时,H=-10dB,满足要求。
其信号流图为:
巴特沃斯高通:
信号流图
巴特沃斯带通:
巴特沃斯带阻
切比雪夫低通
切比雪夫高通
切比雪夫带通
切比雪夫带阻
6总结
通过本次试验,我熟悉了数字滤波的基本概念、数字滤波器的主要技术指标及其物理意义,掌握了巴特沃斯和切比雪夫模拟低通滤波器的设计方法和IIR数字低通滤波器的脉冲响应不变设计法、双线性变换法设计方法,了解了模拟和数字滤波器的频率变换、IIR数字滤波器的直接(优化)设计方法。
7参考资料
史林,赵树杰. 数字信号处理. 北京:科学出版社,2007。