数字信号处理实验 3

上机频谱分析过程及结果图 上机实验三：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 )]使用双线性变换法时，克服了多值映射的关系，避免了频率响应的混叠现象；在零频率附近，频率关系接近于线性关系，高频处有较大的非线性失真。

数字信号处理实验报告姓名：班级：通信学号：实验名称：频域抽样定理验证实验类型：验证试验指导教师:实习日期：2013.频域采样定理验证实验一． 实验目的：1. 加深对离散序列频域抽样定理的理解2.了解由频谱通过IFFT 计算连续时间信号的方法3.掌握用MATLAB 语言进行频域抽样与恢复时程序的编写方法 4、用MATLAB 语言将X(k)恢复为X(z)及X(e jw )。

二． 实验原理：1、1、频域采样定理： 如果序列x(n)的长度为M ，频域抽样点数为N ，则只有当频域采样点数N ≥M 时，才有x N (n)=IDFT[X(k)]=x(n),即可由频域采样X(k)无失真的恢复原序列 x(n)。

2、用X(k)表示X(z)的内插公式：∑-=-----=10111)(1)(N k kNNzWz k X Nz X内插函数： zWzkNNN z 1k111)(-----=ϕ频域内插公式：∑-=-=10)2()()(N K j k Nk X e X πωϕω频域内插函数：e N j N N )21()2sin()2sin(1)(--=ωωωωϕ三． 实验任务与步骤：实验一：长度为26的三角形序列x(n)如图(b)所示，编写MATLAB 程序验证频域抽样定理。

实验二：已知一个时间序列的频谱为X(e jw )=2+4e -jw +6e -j2w +4e -j3w +2e -j4w分别取频域抽样点数N为3、5和10，用IPPT计算并求出其时间序列x(n)，用图形显示各时间序列。

由此讨论原时域信号不失真地由频域抽样恢复的条件。

实验三：由X32(k)恢复X(z)和X(e jw)。

四．实验结论与分析：实验一：源程序：M=26;N=32;n=0:M; %产生M长三角波序列x(n)xa=0:floor(M/2);xb= ceil(M/2)-1:-1:0; xn=[xa,xb];Xk=fft(xn,512); %1024点FFT[x(n)], 用于近似序列x(n)的TFX32k=fft(xn,32); %32点FFT[x(n)]x32n=ifft(X32k); %32点IFFT[X32(k)]得到x32(n)X16k=X32k(1:2:N); %隔点抽取X32k得到X16(K)x16n=ifft(X16k,N/2); %16点IFFT[X16(k)]得到x16(n)subplot(3,2,2);stem(n,xn,'.');box ontitle('(b) 三角波序列x(n)');xlabel('n');ylabel('x(n)');axis([0,32,0,20])k=0:511;wk=2*k/512;subplot(3,2,1);plot(wk,abs(Xk));title('(a)FT[x(n)]');xlabel('\omega/\pi');ylabel('|X(e^j^\omega)|');axis([0,1,0,200])k=0:N/2-1;subplot(3,2,3);stem(k,abs(X16k),'.');box ontitle('(c) 16点频域');xlabel('k');ylabel('|X_1_6(k)|');axis([0,8,0,200])n1=0:N/2-1;subplot(3,2,4);stem(n1,x16n,'.');box ontitle('(d) 16点IDFT[X_1_6(k)]');xlabel('n');ylabel('x_1_6(n)');axis([0,32,0,20])k=0:N-1;subplot(3,2,5);stem(k,abs(X32k),'.');box ontitle('(e) 32点频域采样');xlabel('k');ylabel('|X_3_2(k)|');axis([0,16,0,200])n1=0:N-1;subplot(3,2,6);stem(n1,x32n,'.');box ontitle('(f) 32点IDFT[X_3_2(k)]');xlabel('n');ylabel('x_3_2(n)');axis([0,32,0,20])结果如下所示：实验一分析：序列x(n)的长度M=26，由图中可以看出，当采样点数N=16<M时，x16(n)确实等于原三角序列x(n)以16为周期的周期延拓序列的主值序列。

1. 假设系统用下面差分方程描述：y(n)=x(n)+ay(n-1)假设a=0.7, 0.8, 0.9 ，分别在三种情况下分析系统的频率特性，并打印幅度特性曲线。

解：B=1;A=[1,-0.7];subplot(3,3,3);zplane(B,A);xlabel('实部Re→');ylabel('虚部Im→');title('y(n)=x(n)-0.7y(n-1)传输函数零、极点分布');grid on [H,w]=freqz(B,A,'whole');subplot(3,3,2);plot(w/pi,abs(H),'linewidth',2);grid on;xlabel('\omega/\pi');ylabel('|H(e^j^\omega)|');title('幅频响应特性');axis([0,2,0,6]);subplot(3,3,1);plot(w/pi,angle(H),'linewidth',2);grid on;axis([-0.1,2.1,-3,3]);xlabel('\omega/\pi');ylabel('\phi(\omega)');title('相频响应特性');B=1;A=[1,-0.8];subplot(3,3,6);zplane(B,A);xlabel('实部Re→');ylabel('虚部Im→');title('y(n)=x(n)-0.8y(n-1)传输函数零、极点分布');grid on [H,w]=freqz(B,A,'whole');subplot(3,3,5);plot(w/pi,abs(H),'linewidth',2);grid on;xlabel('\omega/\pi');ylabel('|H(e^j^\omega)|');title('幅频响应特性');axis([0,2,0,6]);subplot(3,3,4);plot(w/pi,angle(H),'linewidth',2);grid on;axis([-0.1,2.1,-3,3]);xlabel('\omega/\pi');ylabel('\phi(\omega)');title('相频响应特性');B=1;A=[1,-0.9];subplot(3,3,9);zplane(B,A);xlabel('实部Re→');ylabel('虚部Im→');title('y(n)=x(n)-0.9y(n-1)传输函数零、极点分布');grid on [H,w]=freqz(B,A,'whole');subplot(3,3,8);plot(w/pi,abs(H),'linewidth',2);grid on;xlabel('\omega/\pi');ylabel('|H(e^j^\omega)|');title('幅频响应特性');axis([0,2,0,6]);subplot(3,3,7);plot(w/pi,angle(H),'linewidth',2);grid on;axis([-0.1,2.1,-3,3]);xlabel('\omega/\pi');ylabel('\phi(\omega)');title('相频响应特性'); 图形如下所示：2.假设系统用下面差分方程描述：y(n) = x(n) +ax(n-1)假设a=0.7, 0.8, 0.9 ，分别在三种情况下分析系统的频率特性，并打印幅度特性曲线。

% Program 5_1% Illustration of DFT Computation% Read in the length N of sequence and the desired % length M of the DFTN = input('Type in the length of the sequence = '); M = input('Type in the length of the DFT = ');% Generate the length-N time-domain sequenceu = [ones(1,N)];% Compute its M-point DFTU = fft(u,M);% Plot the time-domain sequence and its DFTt = 0:1:N-1;stem(t,u)title('Original time-domain sequence')xlabel('Time index n'); ylabel('Amplitude')pausesubplot(2,1,1)k = 0:1:M-1;stem(k,abs(U))title('Magnitude of the DFT samples')xlabel('Frequency index k'); ylabel('Magnitude') subplot(2,1,2)stem(k,angle(U))title('Phase of the DFT samples')xlabel('Frequency index k'); ylabel('Phase')% Program 5_1% Illustration of DFT Computation% Read in the length N of sequence and the desired% length M of the DFTN = input('Type in the length of the sequence = ');M = input('Type in the length of the DFT = ');% Generate the length-N time-domain sequenceu=[ones(1,N) zeros(1,M-N)];% Compute its M-point DFTD=dftmtx(M);U=D*u';% Plot the time-domain sequence and its DFTt = 0:1:M-1;stem(t,u)title('Original time-domain sequence')xlabel('Time index n'); ylabel('Amplitude') pausesubplot(2,1,1)k = 0:1:M-1;stem(k,abs(U))title('Magnitude of the DFT samples')xlabel('Frequency index k'); ylabel('Magnitude') subplot(2,1,2)stem(k,angle(U))title('Phase of the DFT samples')xlabel('Frequency index k'); ylabel('Phase')% Program 5_2% Illustration of IDFT Computation% Read in the length K of the DFT and the desired % length N of the IDFTK = input('Type in the length of the DFT = ');N = input('Type in the length of the IDFT = ');% Generate the length-K DFT sequencek = 0:K-1;V = k/K;% Compute its N-point IDFTv = ifft(V,N);% Plot the DFT and its IDFTstem(k,V)xlabel('Frequency index k'); ylabel('Amplitude') title('Original DFT samples')pausesubplot(2,1,1)n = 0:N-1;stem(n,real(v))title('Real part of the time-domain samples') xlabel('Time index n'); ylabel('Amplitude') subplot(2,1,2)stem(n,imag(v))title('Imaginary part of the time-domain samples') xlabel('Time index n'); ylabel('Amplitude')K = input('Type in the length of the DFT = ');N = input('Type in the length of the IDFT = ');k = 0:K-1;V = k/K;p=[V zeros(1,N-K)];v = conj(dftmtx(N))/N*p';v=v';stem(k,V)xlabel('Frequency index k'); ylabel('Amplitude')title('Original DFT samples')pausesubplot(2,1,1)n = 0:N-1;stem(n,real(v))title('Real part of the time-domain samples')xlabel('Time index n'); ylabel('Amplitude')subplot(2,1,2)stem(n,imag(v))title('Imaginary part of the time-domain samples')xlabel('Time index n'); ylabel('Amplitude')% Program 5_3% Numerical Computation of Fourier transform Using DFTk = 0:15; w = 0:511;x = cos(2*pi*k*3/16);% Generate the length-16 sinusoidal sequence X = fft(x); % Compute its 16-point DFTXE = fft(x,512); % Compute its 512-point DFT% Plot the frequency response and the 16-point DFT samplesplot(k/16,abs(X),'o', w/512,abs(XE))xlabel('\omega/\pi'); ylabel('Magnitude')k=0:15;w=0:511;x=cos(2*pi*k*3/16);X=freqz(x,1,16,'whole');XE=freqz(x,1,512,'whole');plot(k/16,abs(X),'o',w/512,abs(XE));xlabel('\omega/\pi');ylabel('Magnitude');x=input('Type in the first sequence = ');h=input('Type in the second sequence = ');L=length(x)+length(h)-1;XE=fft(x,L);HE=fft(h,L);y1=ifft(XE.*HE);y2=conv(x,h);y3=abs(y1-y2);n=0:L-1;subplot(3,1,1);stem(n,y1);xlabel('Time index n');ylabel('Amplitude');title('Result of DFT-based linear convolution'); subplot(3,1,2);stem(n,y2);xlabel('Time index n');ylabel('Amplitude');title('Result of DFT-based linear convolution'); subplot(3,1,3);stem(n,y3);xlabel('Time index n');ylabel('Amplitude');title('Error sequence');Type in the first sequence = [1 2 0 1]Type in the second sequence = [2 2 1 1]。

实验3 离散时间系统的频域分析一、实验目的（1）了解DFS 、DFT 与DTFT 的联系；加深对FFT 基本理论的理解；掌握用MATLB 语言进行傅里叶变换时常用的子函数；（2）了解离散系统的零极点与系统因果性和稳定性的关系；加深对离散系统的频率响应特性基本概念的理解；熟悉MATLAB 中进行离散系统零极点分析的常用子函数；掌握离散系统幅频响应和相频响应的求解方法。

二、实验内容1. 已知离散时间系统函数为 用matlab 中的函数()432143213.07.05.11.112.01.03.01.02.0--------+-+-++++=zz z z z z z z z H 求该系统的零极点及零极点分布图，并判断系统的因果稳定性。

方法一：利用tf3zp 函数b=[0.2 0.1 0.3 0.1 0.2]; a=[1 -1.1 1.5 -0.7 0.3]; [z,p,k]=tf2zp(b,a); c1=abs(z);c2=angle(z); c3=abs(p);c4=angle(p); polar(c4,c3,'bx') hold onpolar(c2,c1,'ro') disp(z) disp(p)disp(abs(z)) disp(abs(p))90270方法二：利用zplaneb=[0.2 0.1 0.3 0.1 0.2];a=[1 -1.1 1.5 -0.7 0.3];z=roots(b);p=roots(a);zplane(b,a)disp(z)disp(p)disp(abs(z))disp(abs(p))-1-0.500.51-1-0.8-0.6-0.4-0.200.20.40.60.81Real PartI m a g i n a r y P a r t由于极点都在单位圆内，故该系统稳定。

若其收敛域为圆外区域，则系统是因果系统。

2. 已知离散时间系统的系统函数为()432143213.07.05.11.112.01.03.01.02.0--------+-+-++++=z z z z z z z z z H求该系统在π~0频率范围内的绝对幅频响应、相频响应。

实验报告课程名称：数字信号处理实验三：窗函数的特性分析班级：通信1403学生姓名：强亚倩学号：1141210319指导教师：范杰清华北电力大学（北京）一、实验目的分析常用窗函数的时域和频域特性，灵活运用窗函数分析信号频谱和设计FIR数字滤波器。

二、实验原理在确定信号谱分析、随机信号功率谱估计以及FIR数字滤波器设计中，窗函数的选择起着重要的作用。

在信号的频谱分析中，截短无穷长的序列会造成频率泄漏，影响频谱分析的精度和质量。

合理选取窗函数的类型，可以改善泄漏现象。

在FIR数字滤波器设计中，截短无穷长的系统单位脉冲序列会造成FIR滤波器幅度特性的波动，且出现过渡带。

三、实验内容1．分析并绘出常用窗函数的时域特性波形（1）矩形窗函数时域波形及频谱①编程②结果：N=51;w=boxcar(N)Y=fft(w,256);subplot(2,1,1);stem([0:N-1],w);xlabel('w');ylabel('y');title('时域波形');subplot(2,1,2);Y0=abs(fftshift(Y));plot([-128:127],Y0)xlabel('W');ylabel('Y0');title('频谱图形');（2）hanning窗函数时域波形及频谱①编程②结果clear all;clc;n=51;w=hanning(n);y0=fft(w,256);subplot(2,1,1);stem([0:n-1],w)xlabel('n');ylabel('w');title('hanning窗时域波形')subplot(2,1,2);Y=abs(fftshift(y0));plot([-128:127],Y);xlabel('w')ylabel('Y')title('hanning频域波形')（3）哈明窗函数时域波形及频谱①编程clear all;clc;n=51;w=hamming(n);y0=fft(w,256);subplot(2,1,1);stem([0:n-1],w)xlabel('n');ylabel('w');title('hamming窗时域波形')subplot(2,1,2);Y=abs(fftshift(y0));plot([-128:127],Y);xlabel('w')ylabel('Y')title('hamming频域波形')②结果（4）blackman窗函数时域波形及频谱①编程clear all;clc;n=51;w=blackman(n);y0=fft(w,256);subplot(2,1,1);stem([0:n-1],w)xlabel('n');ylabel('w');title('blackman窗时域波形')subplot(2,1,2);Y=abs(fftshift(y0));plot([-128:127],Y);xlabel('w')ylabel('Y')title('blackman频域波形') ②结果（5）battlett窗函数时域波形及频域特性①编程②结果clear all;clc;n=51;w=bartlett(n);y0=fft(w,256);subplot(2,1,1);stem([0:n-1],w)xlabel('n');ylabel('w');title('bartlett窗时域波形')subplot(2,1,2);Y=abs(fftshift(y0));plot([-128:127],Y);xlabel('w')ylabel('Y')title('bartlett频域波形')（6）Kaiser窗函数时域及频域波形①编程clear all;clc;n=51;w=kaiser(n);y0=fft(w,256);subplot(2,1,1);stem([0:n-1],w)xlabel('n');ylabel('w');title('Kaiser时域波形')subplot(2,1,2);Y=abs(fftshift(y0));plot([-128:127],Y);xlabel('w')ylabel('Y')title('Kaiser频域波形')②结果3. 研究凯塞窗(Kaiser)的参数选择对其时域和频域的影响。

数字信号处理高西全实验报告三选择FFT的变换区间N为8和16 两种情况进行频谱分析^p 。

分别打印其幅频特性曲线。

并进行对比、分析^p 和讨论。

（2）对以下周期序列进行谱分析^p 。

选择FFT的变换区间N为8和16 两种情况分别对以上序列进行频谱分析^p 。

分别打印其幅频特性曲线。

并进行对比、分析^p 和讨论。

（3）对模拟周期信号进行谱分析^p选择采样频率，变换区间N=16,32,64 三种情况进行谱分析^p 。

分别打印其幅频特性，并进行分析^p 和讨论。

四、程序码与运行结果(1) 实验程序：1n=[ones(1,4)];M=8;a=1:(M/2); b=(M/2):-1:1; 2n=[a,b];3n=[b,a];1k8=fft(1n,8);1k16=fft(1n,16);2k8=fft(2n,8);2k16=fft(2n,16);3k8=fft(3n,8);3k16=fft(3n,16);以下绘制幅频特性曲线n=0:length(1k8)-1;subplot(3,2,1);stem(n,abs(1k8),#;.#;);label({#;ω/π#;;#;8点DFT[1(n)]#;});ylabel(#;幅度#;);n=0:length(1k16)-1;subplot(3,2,2);stem(n,abs(1k16),#;.#;);label({#;ω/π#;;#;16点DFT[1(n)]#;});ylabel(#;幅度#;); n=0:length(2k8)-1;subplot(3,2,3);stem(n,abs(2k8),#;.#;);label({#;ω/π#;;#; 8点DFT[2(n)]#;});ylabel(#;幅度#;); n=0:length(2k16)-1;subplot(3,2,4);stem(n,abs(2k16),#;.#;);label({#;ω/π#;;#;16点DFT[2(n)]#;});ylabel(#;幅度#;); n=0:length(3k8)-1;subplot(3,2,5);stem(n,abs(3k8),#;.#;);l abel({#;ω/π#;;#; 8点DFT[3(n)]#;});ylabel(#;幅度#;); n=0:length(3k16)-1;subplot(3,2,6);stem(n,abs(3k16),#;.#;);label({#;ω/π#;;#;16点DFT[3(n)]#;});ylabel(#;幅度#;); 图形：（2）实验程序：n=0:7;4n=cos(pi/4n);4k8=fft(4n,8);subplot(2,2,1);stem(2n/8,abs(4k8),#;.#;);label({#;ω/π#;;#;8点DFT[4(n)]#;});ylabel(#;幅度#;); 5n=cos(pi/4n)+cos(pi/8n);5k8=fft(5n,8);subplot(2,2,2);stem(2n/8,abs(5k8),#;.#;);label({#;ω/π#;;#;8点DFT[5(n)]#;});ylabel(#;幅度#;); n=0:15;4n=cos(pi/4n);5n=cos(pi/4n)+cos(pi/8n);4k16=fft(4n,16);subplot(2,2,3);stem(2n/16,abs(4k16),#;.#;);label({#;ω/π#;;#;16点DFT[4(n)]#;});ylabel(#;幅度#;); 5k16=fft(5n,16);subplot(2,2,4);stem(2n/16,abs(5k16),#;.#;);label({#;ω/π#;;#;16点DFT[5(n)]#;});ylabel(#;幅度#;); 图形：（3）实验代码：Fs=64;T=1/Fs;N=16;n=0:N-1;6nT=cos(8pinT)+cos(16pinT)+cos(20pinT);6k16=fft(6nT);6k16=fftshift(6k16);Tp=NT;F=1/Tp;k=-N/2:N/2-1;fk=kF;subplot(3,1,1);stem(fk,abs(6k16),#;.#;);label({#;f(Hz)#;;#;16点DFT[6(nT)]#;});ylabel(#;幅度#;); N=32;n=0:N-1;6nT=cos(8pinT)+cos(16pinT)+cos(20pinT);6k32=fft(6nT,32);6k32=fftshift(6k32);Tp=NT;F=1/Tp;k=-N/2:N/2-1;fk=kF;subplot(3,1,2);stem(fk,abs(6k32),#;.#;);label({#;f(Hz)#;;#;32点DFT[6(nT)]#;});ylabel(#;幅度#;); N=64;n=0:N-1;6nT=cos(8pinT)+cos(16pinT)+cos(20pinT);6k64=fft(6nT,64);6k64=fftshift(6k64);Tp=NT;F=1/Tp;k=-N/2:N/2-1;fk=kF;subplot(3,1,3);stem(fk,abs(6k64),#;.#;);label({#;f(Hz)#;;#;64点DFT[6(nT)]#;});ylabel(#;幅度#;);图形：五、实验总结1．结论用DFT对信号进行谱分析^p 时，重点关注频谱分辨率和分析^p 误差，频谱分辨率F=1/Tp=Fs/N，可以依据此等式来选择FFT的变换区间N，而误差主要来自于用FFT作频谱分析^p 时，得到的是离散谱，而当信号是非周期信号时，应该得到连续谱，只有当N较大时，用FFT做出来的离散谱才接近于连续谱，因此N要适当选择大一些。

数字信号处理实验三离散时间傅里叶变换DTFT及IDTFT一、实验目的：(1)通过本实验,加深对DTFT和IDFT的理解;(2)熟悉应用DTFT对典型信号进行频谱分析的方法;(3)掌握用MATLAB进行离散时间傅里叶变换及其逆变换的方法;二、实验内容：1自己生成正弦序列如矩形序列,正弦序列,指数序列等,对其进行频谱分析,观察其时域波形和频域的幅频特性;记录实验中观察到的现象,绘出相应的时域序列和幅频特性曲线;矩形序列：程序：M=10;N=2M+1;T=;n=-4M:4M;x=zeros1,3M,ones1,N,zeros1,3M;w=-15::15+1e-10;X=sinNwT./sinwT;subplot1,3,1;stemn,x,'.';axis-20,20,,,grid onxlabel'n',title'a序列幅度'subplot1,3,2,plotw,X,grid onxlabel'\Omega',title'b幅频特性'subplot1,3,3,plotw,X,grid onv=axis;axis-pi/T,pi/T,v3,v4;xlabel'\Omega',title'c横轴放大后幅频特性' setgcf,'color','w'正弦序列：程序:n=-10:10; x=sinnpi;k=-200:200; w=pi/100k;X=xexp-jpi/100.^n'k; magX=absX;angX=angleX;subplot3,1,1;stemn,x,'.k';title'xn=sinπn';subplot3,1,2;plotw/pi,magX,'.k';title'Xe^jw幅度谱';subplot3,1,3;plotw/pi,angX,'.k';title'Xe^jw相位谱';n=-10:10; x=sinnpi;k=-200:200;w=pi/100k;X=xexp-jpi/100.^n'k;magX=absX;angX=angleX;subplot3,1,1;stemn,x,'.k';title'xn=sinπn';subplot3,1,2;plotw/pi,magX,'.k'; title'Xe^jw幅度谱'; subplot3,1,3;plotw/pi,angX,'.k'; title'Xe^jw相位谱';波形如下:指数序列：程序：n=-5:5;x=.^n;k=-200:200;w=pi/100k;X=xexp-jpi/100.^n'k;magX=absX;angX=angleX;subplot2,1,1;plotw/pi,magX,'k';grid;axis-2,2,0,15xlabel'frequency in units of\pi';ylabel'｜x｜'gtext'Magnitde Part'subplot2,1,2;plotw/pi,angX,'k'/pi,grid;axis-2,2,-4,4xlabel'frequency in units of\pi';ylabel'radians\pi' gtext'Angle Part';2.对于理想的低通,高通滤波器,用IDTFT 求出它的逆变换所对应得离散时间序列;记录实验中观察到的现象,绘出相应的时域序列曲线;要求滤波器的截至频率可由用户在MATLAB 界面自行输入;程序：wc=pi;n=-10:10+1e-10;hd=sinnwc./npi;subplot1,2,1;plot-pi,-wc,-wc,wc,wc,pi,0,0,1,1,0,0xlabel'频率1/秒';ylabel'幅度';axis-pi,pi,,,grid onsubplot1,2,2;stemn,hd,grid onxlabel'n';ylabel'序列';axis-10,10,wc,wcsetgcf,'color','w'三、思考题离散时间信号的频谱分辨率在实验中能体现出来吗实序列的DTFT具有对称性吗若是,如何体现出来答：能,实序列的DTFT具有对称性;离散时间信号的频谱中,频谱分辨率体现在相同的坐标系下面,能表现信号的范围,当表现的范围越大,其分辨率越高。

实验三：用FFT 对信号做频谱分析1 实验目的学习用FFT 对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便正确应用FFT 。

2 实验原理用FFT 对信号做频谱分析是学习数字信号处理的重要内容。

经常需要进行谱分析的信号是模拟信号和时域离散信号。

对信号进行谱分析的重要问题是频谱分辨率D 和分析误差。

频谱分辨率直接和FFT 的变换区间N 有关，因为FFT 能够实现的频率分辨率是N /2π，因此要求D N ≤/2π。

可以根据此式选择FFT 的变换区间N 。

误差主要来自于用FFT 作频谱分析时，得到的是离散谱，而信号（周期信号除外）是连续谱，只有当N 较大时离散谱的包络才能逼近于连续谱，因此N 要适当选择大一些。

周期信号的频谱是离散谱，只有用整数倍周期的长度作FFT ，得到的离散谱才能代表周期信号的频谱。

如果不知道信号周期，可以尽量选择信号的观察时间长一些。

对模拟信号进行谱分析时，首先要按照采样定理将其变成时域离散信号。

如果是模拟周期信号，也应该选取整数倍周期的长度，经过采样后形成周期序列，按照周期序列的谱分析进行。

3 实验步骤及内容（1）对以下序列进行谱分析。

⎪⎩⎪⎨⎧≤≤-≤≤-=⎪⎩⎪⎨⎧≤≤-≤≤+==其它nn n n n n x 其它n n n n n n x n R n x ,074,330,4)(,074,830,1)()()(3241 选择FFT 的变换区间N 为8和16 两种情况进行频谱分析。

分别打印其幅频特性曲线。

并进行对比、分析和讨论。

（2）对以下周期序列进行谱分析。

4()c o s 4x n n π= 5()c o s (/4)c o s (/8)x n n n ππ=+ 选择FFT 的变换区间N 为8和16 两种情况分别对以上序列进行频谱分析。

分别打印其幅频特性曲线。

并进行对比、分析和讨论。

（3）对模拟周期信号进行谱分析6()cos8cos16cos20x t t t t πππ=++选择 采样频率Hz F s 64 ，变换区间N=16,32,64 三种情况进行谱分析。

实验三 抽样一、实验原理抽样过程基于两个基本原理：混叠合重建，涉及正弦波合现行调频信号的混叠。

二、 实验内容.3.3.1抽样引起的混叠 由于在MATLAB 中不能产生模拟信号，实验需要做实时t 轴的仿真。

因此，把仿真时的△t 与所研究的抽样周期Ts 明确地区分开始很重要的。

1. 正弦信号混叠对连续时间正弦信号考虑下面表达式：()()φπ+=t f t 02sin x可以按抽样频率ST s f 1=对()t x 抽样来获得离散时间信号[]()()⎪⎪⎭⎫ ⎝⎛+=====φπn f f t x t x s f n t nT t ss 0/2sin n x 2. 实验内容以不同组合的0f 和s f 绘出[]n x ，可以说明混叠问题。

以下，取抽样频率s f =8kHz.a. 首先，绘出一个被抽样的正弦波的单图。

令正弦波的频率为300Hz ，然后在10ms 长间隔上抽样。

相位φ可以任意指定。

使用stem 绘出产生的离散时间信号。

因为是用眼睛实现重建可视化信号包络，应该很容易看到正弦信号的轮廓。

b. 如果必要，使用plot 绘图。

在这种情况下，点用直线段连接 起来， 这样正弦信号的特点应该是明显的。

用直线段连接信号 样本是自离散时间样本产生连续时间信号的“信号重建”的 一 种方式。

它 不是抽样定理所说的理想重建，但对大多数的情形它已经是足够好，很有用。

c. 把正弦的频率从100Hz 变至475Hz,每次增加125Hz ，如在()a 部 分中那样，绘出一系列相应的图。

注意，正如所预期的那样， 显 现的正弦信号的频率在逐渐增加。

最好用subplot 指令把四 个图放在同一屏上。

d. 把正弦的频率从7525Hz 变至7900Hz ，每次增加125Hz ，正 如在()c 中那样，另外绘出一系列相应的图。

注意现在显现的正弦信号的频率在逐渐减少。

解释这一现象。

e.把正弦的频率从32100Hz变至32475Hz，每次增加125Hz，再次绘出一些列类似的图。

一、实验目的1．了解工程上两种最常用的变换方法：脉冲响应不变法和双线性变换法。

2．掌握双线性变换法设计IIR 滤波器的原理及具体设计方法，熟悉用双线性设计法设计低通、带通和高通IIR 数字滤波器的计算机程序.3．观察用双线性变换法设计的滤波器的频域特性，并与脉冲响应不变法相比较,了解双线性变换法的特点。

4．熟悉用双线性变换法设计数字Butterworth 和Chebyshev 滤波器的全过程。

5．了解多项式乘积和多项式乘方运算的计算机编程方法。

二、实验原理与方法从模拟滤波器设计IIR 数字滤波器具有四种方法：微分－差分变换法、脉冲响应不变法、双线性变换法、z 平面变换法。

工程上常用的是其中的两种:脉冲响应不变法、双线性变换法。

脉冲响应不变法需要经历如下基本步骤：由已知系统传输函数H(S)计算系统冲激响应h(t)；对h(t)等间隔采样得到h （n ）=h （n T)；由h （n ）获得数字滤波器的系统响应H （Z)。

这种方法非常直观，其算法宗旨是保证所设计的IIR 滤波器的脉冲响应和模拟滤波器的脉冲响应在采样点上完全一致。

而双线性变换法的设计准则是使数字滤波器的频率响应与参考模拟滤波器的频率响应相似。

脉冲响应不变法一个重要的特点是频率坐标的变换是线性的（），其确定是有频谱的周期延拓效应，存在频谱混叠的现象。

为了克服脉冲响应不变法可能产生的频谱混叠，提出了双线性变换法，它依靠双线性变换式：， ， 其中 ,建立其S 平面和Z 平面的单值映射关系，数字域频率和模拟域频率的关系是： , （3－1） 由上面的关系式可知,当时，终止在折叠频率处,整个轴单值的对应于单位圆的一周。

因此双线性变换法不同于脉冲响应不变法，不存在频谱混叠的问题。

从式（3－1）还可以看出，两者的频率不是线性关系。

这种非线性关系使得通带截至频率、过渡带的边缘频率的相对位置都发生了非线性畸变。

这种频率的畸变可以通过预畸变来校正。

用双线性变换法设计数字滤波器时，一般总是先将数字滤波器的个临界频率经过式（3－1）的频率预畸变,求得相应参考模拟滤波器的个临界频率，然后设计参考模拟滤波器的传递函数，最后通过双T Ω=ω1111--+-=z z s s s z -+=11Ω+=j s σωj re z =)2/(ωtg =Ω)(2Ω=arctg ω∞→Ωωπω=Ωj线性变换式求得数字滤波器的传递函数。

数字信号处理实验三数字信号处理实验三是针对数字信号处理课程的一项实践性任务。

本实验旨在通过实际操作，加深对数字信号处理理论的理解，并培养学生的实验能力和问题解决能力。

在本实验中，我们将学习和实践以下内容：1. 实验目的本实验的目的是通过使用MATLAB软件进行数字信号处理，加深对数字信号处理基本概念和算法的理解，掌握数字信号的采样、量化、滤波等基本操作。

2. 实验器材在本实验中，我们将使用以下器材：- 个人计算机- MATLAB软件3. 实验步骤本实验的具体步骤如下：步骤一：信号生成首先，我们需要生成一个模拟信号，可以是正弦信号、方波信号或其他类型的信号。

在MATLAB中，我们可以使用相关函数生成这些信号。

生成信号的目的是为了后续的数字信号处理操作提供输入。

步骤二：信号采样在本步骤中，我们将对生成的模拟信号进行采样。

采样是指在一定的时间间隔内对信号进行离散化处理，得到离散时间上的信号序列。

在MATLAB中，我们可以使用采样函数对信号进行采样。

步骤三：信号量化在本步骤中，我们将对采样后的信号进行量化。

量化是指将连续的信号离散化为一组离散的幅值。

在MATLAB中，我们可以使用量化函数对信号进行量化。

步骤四：信号滤波在本步骤中，我们将对量化后的信号进行滤波。

滤波是指通过一系列滤波器对信号进行处理，以去除不需要的频率成分或噪声。

在MATLAB中，我们可以使用滤波函数对信号进行滤波。

步骤五：信号重构在本步骤中，我们将对滤波后的信号进行重构。

重构是指将离散化的信号恢复为连续的信号。

在MATLAB中，我们可以使用重构函数对信号进行重构。

步骤六：信号分析在本步骤中，我们将对重构后的信号进行分析。

分析是指对信号的频谱、功率等特性进行分析，以了解信号的特点和性能。

在MATLAB中，我们可以使用分析函数对信号进行分析。

4. 实验结果在完成以上步骤后，我们可以得到经过数字信号处理的结果。

这些结果可以是经过采样、量化、滤波和重构后的信号波形，也可以是信号的频谱、功率等特性。

实验三快速Fourier变换(FFT)及其应用一、实验目的1．在理论学习的基础上，通过本实验，加深对FFT的理解，熟悉MATLAB 中的有关函数。

2．应用FFT对典型信号进行频谱分析。

3. 了解应用FFT进行信号频谱分析过程中可能出现的问题，以便在实际中正确应用FFT。

4．应用FFT实现两个序列的线性卷积和方法。

二、实验原理与方法上机实验内容：(1)、观察高斯序列的时域和幅频特性，固定信号xa(n)中参数p=8，改变q的值，使q分别等于2，4，8，观察它们的时域和幅频特性，了解当q取不同值时，对信号序列的时域幅频特性的影响；固定q=8，改变p，使p分别等于8，13，14，观察参数p变化对信号序列的时域及幅频特性的影响，观察p等于多少时，会发生明显的泄漏现象，混叠是否也随之出现？记录实验中观察到的现象，绘出相应的时域序列和幅频特性曲线。

程序：function [x, F]=gauss(p,q);n=0:15;x(n+1)=exp(-(n+1-p).^2/q);F=fft(x);endclear all;figure(1)[x1,F1]=gauss(8,2);n=0:15;subplot(3,2,1);plot(n,x1);text(6,0.2,'p=8,q=2');grid on;xlabel('n');ylabel('时域');subplot(3,2,2);plot(abs(F1));text(7.5,2,'p=8,q=2');grid on;xlabel('k');ylabel('频域');[x2,F2]=gauss(8,4);n=0:15;subplot(3,2,3);plot(n,x2);text(6,0.2,'p=8,q=4');grid on;xlabel('n');ylabel('时域');subplot(3,2,4);plot(abs(F2));text(7.5,2,'p=8,q=4'); grid on;xlabel('k');ylabel('频域');[x3,F3]=gauss(8,8);n=0:15;subplot(3,2,5);plot(n,x3);text(6,0.2,'p=8,q=8'); grid on;xlabel('n');ylabel('时域');subplot(3,2,6);plot(abs(F3));text(7.5,2,'p=8,q=8'); grid on;xlabel('k');ylabel('频域');figure(2)[x4,F4]=gauss(8,8);n=0:15;subplot(3,2,1);plot(n,x4);text(6,0.2,'p=8,q=8'); grid on;xlabel('n');ylabel('时域');subplot(3,2,2);plot(abs(F4));text(7.5,2,'p=8,q=8'); grid on;xlabel('k');ylabel('频域');[x5,F5]=gauss(13,8);n=0:15;subplot(3,2,3);plot(n,x5);text(6,0.2,'p=13,q=8'); grid on;xlabel('n');ylabel('时域');subplot(3,2,4);plot(abs(F5));text(7.5,2,'p=13,q=8'); grid on ;xlabel('k');ylabel('频域');[x6,F6]=gauss(14,8); n=0:15;subplot(3,2,5);plot(n,x6);text(6,0.2,'p=14,q=8'); grid on ;xlabel('n');ylabel('时域'); subplot(3,2,6);plot(abs(F6));text(7.5,2,'p=14,q=8'); grid on ;xlabel('k'); ylabel('频域');5101500.51p=8,q=2n时域0510152024p=8,q=2k频域5101500.51p=8,q=4n时域0510152024p=8,q=4k频域5101500.51p=8,q=8n时域51015200510p=8,q=8k频域5101500.51p=8,q=8n时域51015200510p=8,q=8k频域5101500.51p=13,q=8n时域510152005p=13,q=8k频域5101500.51p=14,q=8n时域510152005p=14,q=8k频域(2)、观察衰减正弦序列xb(n)的时域和幅频特性，a=0.1，f=0.0625，检查谱峰出现位置是否正确，注意频谱的形状，绘出幅频特性曲线，改变f ，使f 分别等于0.4375和0.5625，观察这两种情况下，频谱的形状和谱峰出现位置，有无混叠和泄漏现象？说明产生现象的原因。

一、实验室名称：数字信号处理实验室 二、实验项目名称：采样的时域及频域分析 三、实验原理：1、采样的概念：采样是将连续信号变化为离散信号的过程。

1. A 、理想采样：即将被采样信号与周期脉冲信号相乘B 、实际采样：将被采样信号与周期门信号相乘，当周期门信号的宽度很小，可近似为周期脉冲串。

根据傅里叶变换性质000()()()()ˆˆ()()()()()()(())FTFTa a T n n FTa a T a T a an n x t X j T j xt x t T x nT t nT X j Xj n ωδωδδδω=+∞=+∞=-∞=-∞←−→Ω←−→Ω==-←−→Ω=Ω-Ω∑∑式中T 代表采样间隔，01TΩ=由上式可知：采样后信号的频谱是原信号频谱以0Ω为周期的搬移叠加 结论：时域离散化，频域周期化；频谱周期化可能造成频谱混迭。

)(t T δ^T ^)tC 、低通采样和Nyquist 采样定理设()()a a x t X j ⇔Ω且()0,2a M M X j f πΩ=Ω>Ω=当，即为带限信号。

则当采样频率满足2/22s M M f f π≥Ω=时,可以从采样后的^()()()a assn x t x nT t nT δ∞=-∞=-∑信号无失真地恢复()ax t 。

称2Mf为奈奎斯特频率，12N M T f =为奈奎斯特间隔。

注意：实际应用中，被采信号的频谱是未知的，可以在ADC 前加一个滤波器（防混迭滤波器）。

2、低通采样中的临界采样、欠采样、过采样的时域及频域变化情况。

低通采样中的临界采样是指在低通采样时采样频率2s M f f = 低通采样中的欠采样是指在低通采样时采样频率2s M f f ≤ 低通采样中的欠采样是指在低通采样时采样频率2s M f f ≥ 设一带限信号的频谱如下：ˆ()a xt )(ˆΩj X a()a G j Ω0 m-ΩΩm Ω0T TT-ΩTΩ（1）临界采样（2）过采样（3）欠采样由上图可知，当为临界采样和过采样时，理论上可以无失真的恢复采样信号，但是实际在临界采样时，由于实际滤波器的性能限制，无法无失真的恢复，在欠采样时只能部分恢复原信号的频谱特性。

实验三 用双线性变换法设计IIR 数字滤波器1. 实验目的(1) 熟悉用双线性变换法设计IIR 数字滤波器的原理与方法。

(2) 掌握数字滤波器的计算机仿真方法。

(3) 通过观察对实际心电图信号的滤波作用， 获得数字滤波的感性知识。

2. 实验内容(1) 用双线性变换法设计一个巴特沃斯低通IIR 数字滤波器。

设计指标参数为：在通带内频率低于0.2π时，最大衰减小于1dB ；在阻带内[0.3π, π] 频率区间上，最小衰减大于15dB 。

(2) 以 0.02π为采样间隔， 打印出数字滤波器在频率区间［0, π/2］上的幅频响应特性曲线。

(3) 用所设计的滤波器对实际心电图信号采样序列(在本实验后面给出)进行仿真滤波处理，并分别打印出滤波前后的心电图信号波形图， 观察总结滤波作用与效果。

3.实验原理为了克服用脉冲响应不变法产生频谱混叠现象，可以采用非线性频率压缩方法（正切变换），从s 平面映射到s1平面，再从s1平面映射到z 平面,即实现了双线性变换。

4. 实验步骤(1) 复习有关巴特沃斯模拟滤波器设计和用双线性变换法设计IIR 数字滤波器的内容， 按照例 6.4.2， 用双线性变换法设计数字滤波器系统函数H(z)。

例 6.4.2 中已求出满足本实验要求的数字滤波器系统函数：(2－1)161212120.0007378(1)()(1 1.2680.705)(1 1.01060.3583)(10.9040.215)z H z zz zz z z -------+=-+-+-+31()k K H z ==∏(2－2)A=0.09036B1=1.2686, C1=－0.7051 B2=1.0106, C2=－0.3583 B3=0.9044, C3=－0.2155由(2－1)式和(2－2)式可见，滤波器H(z)由三个二阶滤波器H1(z)，H2(z)和H3(z)级联组成，如图 2－1 所示。

(2) 编写滤波器仿真程序，计算H(z)对心电图信号采样序列x(n)的响应序列y(n)。

数字信号处理实验三报告数字信号处理上机实验报告实验一系统响应及系统稳定性一、实验目的（1）掌握求系统响应的方法。

（2）掌握时域离散系统的时域特性。

（3）分析、观察及检验系统的稳定性。

二、实验内容（1）给定一个低通滤波器的差分方程为y(n)=0.05x(n)+0.05x(n-1)+0.9y(n-1)输入信号x1(n)=R8(n)x2(n)=u(n)(a) 分别求出系统对x1(n)=R8(n) 和x2(n)=u(n)的响应序列，并画出其波形。

(b) 求出系统的单位冲响应，画出其波形。

实验程序：A=[1,-0.9];B=[0.05,0.05]; %%系统差分方程系数向量 B 和 Ax1n=[1 1 1 1 1 1 1 1 zeros(1,50)]; %产生信号 x1(n)=R8(n)x2n=ones(1,8); %产生信号 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)+2.5δ(n-1)+2.5δ(n-2)+δ(n-3)用线性卷积法分别求系统h1(n)和h2(n)对x1(n)=R8(n)的输出响应，波形。

实验三：离散时间信号的频域分析一．实验目的1．在学习了离散时间信号的时域分析的基础上，对这些信号在频域上进行分析，从而进一步研究它们的性质。

2．熟悉离散时间序列的3种表示方法：离散时间傅立叶变换（DTFT）,离散傅立叶变换（DFT）和Z变换。

二．实验相关知识准备1．用到的MATLAB命令运算符和特殊字符：< > .* ^ .^语言构造与调试：error function pause基本函数：angle conj rem数据分析和傅立叶变换函数：fft ifft max min工具箱：freqz impz residuez zplane三．实验内容1．离散傅立叶变换在MATLAB中，使用fft可以很容易地计算有限长序列x[n]的离散傅立叶变换。

此函数有两种形式：y=fft(x)y=fft(x,n) 求出时域信号x的离散傅立叶变换n为规定的点数，n的默认值为所给x的长度。

当n取2的整数幂时变换的速度最快。

通常取大于又最靠近x的幂次。

（即一般在使用fft函数前用n=2^nextpow2(length(x))得到最合适的n)。

当x的长度小于n时，fft函数在x的尾部补0，以构成长为n点数据。

当x的长度大于n时，fft函数将序列x截断，取前n点。

一般情况下，fft求出的函数多为复数，可用abs及angle分别求其幅度和相位。

注意：栅栏效应，截断效应（频谱泄露和谱间干扰），混叠失真例3－1：fft函数最通常的应用是计算信号的频谱。

考虑一个由100hz和200hz正弦信号构成的信号，受零均值随机信号的干扰，数据采样频率为1000hz。

通过fft函数来分析其信号频率成分。

t=0:0.001:1;%采样周期为0.001s，即采样频率为1000hzx=sin(2*pi*100*t)+sin(2*pi*200*t)+1.5*rand(1,length(t));%产生受噪声污染的正弦波信号subplot(2,1,1);plot(x(1:50));%画出时域内的信号y=fft(x,512);%对x进行512点的fftf=1000*(0:256)/512;%设置频率轴（横轴）坐标，1000为采样频率subplot(2,1,2);plot(f,y(1:257));%画出频域内的信号实验内容3－2：频谱泄漏和谱间干扰假设现有含有三种频率成分的信号x(t)=cos(200πt)+sin(100πt)+cos(50πt)用DFT分析x(t)的频谱结构。

数字信号处理实验三复指数序列的绘图一、实验目的:1.充分熟悉复指数函数exp的使用；2.熟悉复指数函数的实部、虚部、振幅、相位的计算；3.能够画出复指数函数实部、虚部、振幅、相位的图形。

二、实验步骤:1.用help查找exp函数的使用情况；2.编辑并生成函数exp.m(单位脉冲序列)n=[0:1:20];alpha=-0.1+0.5j;x=exp(alpha*n);subplot(2,2,1);stem(n,real(x));title('实部');xlabel('n')subplot(2,2,3);stem(n,imag(x));title('虚部');xlabel('n')subplot(2,2,2);stem(n,abs(x));title('振幅');xlabel('n')subplot(2,2,4);stem(n,(180/pi)*angle(x));title('相位');xlabel('n')3.先运行exp.m, 画出图形，分析图形的结果。

分析：当n在0到20的范围递增时，函数的振幅呈递减趋势。

当n在n+12k到n+6+12k的范围递增时，函数的相位呈递增趋势。

4.当alpha=-0.1+0.5j 时，计算x的共轭对称分量和共轭反对称分量.也分实部、虚部、振幅、相位分别画图。

5. 再改变alpha=1.5, alpha=-2, 画出图形，分析图形的结果。

三、实验报告要求：1. 实验目的2. 实验步骤3.14个图形的结果，并得出结论或体会。

4.思考题：试用公式表示x=exp(alpha*n)的实部、虚部、振幅、相位。