《数字信号处理》上机全部源代码调试通过,完整版

合集下载

课程名称数字信号处理教程-课程代码0817

天津市高等教育自学考试课程考试大纲课程名称: 数字信号处理课程代码:0817、4701、4659第一部分课程性质与设置目标一、课程性质与特点数字信号处理是高等教育自学考试电子信息工程专业所开设的专业基础课之一,它是一门理论性较强的课程。

本课程系统地讨论了数字信号处理的基本理论、基本分析方法和基本实现方法。

深入浅出，条理清楚，具有实例分析。

二、课程设置的目的和要求设置本课程，为了使学生牢固掌握数字信号处理的基本理论、基本分析方法和基本实现方法。

通过本课程的学习，要求学生掌握数字信号处理的基本理论，训练用基本分析方法和基本实现方法分析和解决实际问题，从而提高学生分析问题和解决问题的能力。

三、与本专业其它课程的关系数字信号处理是电子信息工程专业学生的专业基础课,它与本专业的其他课程有着密切的关系。

电路理论、信号与系统、积分变换是本课程的基础。

第二部分课程内容与考核目标第一章离散时间信号与系统一、学习的目的和要求通过本章的学习，掌握连续信号的抽样、离散时间信号的概念和线性移不变系统的概念以及描述该系统的数学模型—常系数线性差分方程。

二、考核知识点与考核目标（一）离散时间信号-序列（重点）理解：离散时间信号-序列的概念、图形；应用：序列的运算，包括：移位、翻褶、和、积、累加、差分、卷积；常用序列、序列的周期性、用单位抽样序列来表示任意序列。

（二）线性移不变系统（重点）应用：线性移不变系统的概念、性质、因果系统、系统稳定性的概念；（三）常系数线性差分方程(一般)应用：序列域求解法；（四）连续时间信号抽样（重点）应用：理想抽样、奈奎斯特抽样定理、实际抽样。

第二章 Z变换一、学习的目的和要求通过本章的学习，复习Z变换的定义与收敛域，基本性质和定理，Z反变换，为在变换域分析序列和系统做准备。

理解序列的Z变换与连续信号的拉普拉斯变换、傅里叶变换的关系，序列的傅里叶变换和一些对称性质，离散系统的系统函数、系统的频率响应。

数字信号处理课件第10章上机实验

Y
改变信号序列否?
N 调用系统单位脉冲响应序列产生子程序，h (求n )
调用傅氏变换数值计算子程序，H求(e jω k )
两次调用绘图子程序，分别绘h制(n )，H (e jω k) 图形
Y
改变h (n )否 ?
2019/10/30
课件
18
2019/10/30
第10章上机实验
x1(n ) R 4(n )
n 1,
x2(n ) 8 n

0
0n3 4n7
其它n
4 n 0 n 3
x3(n ) n 3 4 n 7

0

x4(n) cos 4 n
x(n)xa(nT)
x(n)的傅里叶变换为

X(ej) x(n)ejn
2019/10/30
n
课件
(10.3.5)
6
第10章上机实验
比较(10.3.5)和(10.3.4)可知
X ^a(j)X(ej)T
(10.3.6)
在数字计算机上观察分析各种序列的频域特性，
通常对X(ejω)在［0, 2π］上进行M点采样来观察分
课件
第10章上机实验
21
第10章上机实验
3. 上机实验内容 (1) 对 2 中所给出的信号逐个进行谱分析。下面给出针对各信号的FFT变换区间N以及对连续信号x6(t)的采样频率fs，供实验时参考。 x1(n), x2(n), x3(n), x4(n), x5(n): N=8, 16=x6(t): fs=64(Hz), N=16, 32, 64 (2) 令x(n)=x4(n)+x5(n)，用FFT计算 8 点和 16 点离散傅里叶变换，

数字信号处理上机实验代码

数字信号处理上机实验代码文件名：tstem.m（实验一、二需要）程序：f unction tstem(xn,yn)%时域序列绘图函数%xn:被绘图的信号数据序列，yn:绘图信号的纵坐标名称（字符串）n=0:length(xn)-1;stem(n,xn,'.');xlabel('n');ylabel('yn');axis([0,n(end),min(xn),1.2*max(xn)]);文件名：tplot.m（实验一、四需要）程序：function tplot(xn,T,yn)%时域序列连续曲线绘图函数%xn:信号数据序列,yn:绘图信号的纵坐标名称（字符串）%T为采样间隔n=0;length(xn)-1;t=n*T;plot(t,xn);xlabel('t/s');ylabel(yn);axis([0,t(end),min(xn),1.2*max(xn)]);文件名：myplot.m（实验一、四需要）%(1)myplot;计算时域离散系统损耗函数并绘制曲线图。

function myplot(B,A)%B为系统函数分子多项式系数向量%A为系统函数分母多项式系数向量[H,W]=freqz(B,A,1000)m=abs(H);plot(W/pi,20*log10(m/max(m)));grid on;xlabel('\omega/\pi');ylabel('幅度(dB)')axis([0,1,-80,5]);title('损耗函数曲线');文件名：mstem.m（实验一、三需要）程序：function mstem(Xk)%mstem(Xk)绘制频域采样序列向量Xk的幅频特性图M=length(Xk);k=0:M-1;wk=2*k/M;%产生M点DFT对应的采样点频率（关于pi归一化值）stem(wk,abs(Xk),'.');box on;%绘制M点DFT的幅频特性图xlabel('w/\pi');ylabel('幅度');axis([0,2,0,1.2*max(abs(Xk))]);文件名：mpplot.m（实验一需要）%(2)mpplot;计算时域离散系统损耗函数和相频特性函数，并绘制曲线图。

《数字信号处理》上机全部源代码调试通过-完整版

《数字信号处理》上机全部源代码调试通过，完整版（高西全，第四版）实验一%实验1：系统响应及系统稳定性close all;clear all%调用fliter解差分方程，由系统对un的响应判断稳定性%内容1：调用filter解差分方程，由系统对u(n)的响应判断稳定性A=[1,-0.9];B=[0.05,0.05];x1n=[1 1 1 1 1 1 1 1 zeros(1,50)];x2n=ones(1,128);hn=impz(B,A,58);subplot(2,2,1);y='h(n)';tstem(hn,y);title('(a) 系统单位脉冲响应h(n)')y1n=filter(B,A,x1n);subplot(2,2,2);y='y1(n)';tstem(y1n,y);title('(b) 系统对R8(n)的响应y1(n)')y2n=filter(B,A,x2n);subplot(2,2,4);y='y2(n)';tstem(y2n, y);title('(c) 系统对u(n)的响应y2(n)')y1n=filter(B,A,x1n);subplot(2,2,2);y='y1(n)';tstem(y1n,y);title('(b) 系统对R8(n)的响应y1(n)')y2n=filter(B,A,x2n);subplot(2,2,4);y='y2(n)';tstem(y2n, y);title('(c) 系统对u(n)的响应y2(n)')%内容2：调用conv函数计算卷积x1n=[1 1 1 1 1 1 1 1]; %产生信号x1n=R8nh1n=[ones(1,10) zeros(1,10)];h2n=[1 2.5 2.5 1 zeros(1,10)]y21n=conv(h1n,x1n);y22n=conv(h2n,x1n);figure(2)subplot(2,2,1);y='h1(n)';tstem(h1n,y);%调用函数tstem绘图title('(d) 系统单位脉冲响应h1(n)')subplot(2,2,2);y='y21(n)';tstem(y21n,y);title('(e) h1(n)与R8(n)的卷积y21(n)')subplot(2, 2,3); y='h2(n)';tstem(h2n,y); %调用函数tstem绘图title('(f) 系统单位脉冲响应h2(n)')subplot(2, 2, 4); y='y22(n)';tstem(y22n, y);title('(g) h2(n)与R8(n)的卷积y22(n)')%=====================================%内容3：谐振器分析un=ones(1, 256); %产生信号unn=0:255;xsin=sin(0.014*n)+sin(0.4*n); %产生正弦信号A=[1,-1.8237,0.9801];B=[1/100.49,0,-1/100.49];%系统差分方程系数向量B和Ay31n=filter(B,A,un); %谐振器对un的响应y31n y32n=filter(B,A,xsin);%谐振器对正弦信号的响应y32nfigure(3)subplot(2,1,1);y='y31(n)';tstem(y31n,y)title('(h) 谐振器对u(n)的响应y31(n)')subplot(2,1,2);y='y32(n)';tstem(y32n,y);title('(i) 谐振器对正弦信号的响应y32(n)')function tstem(xn,yn)n = 0:length(xn)-1;stem(n,xn,'.');xlabel('n');ylabel('yn');%xlabel('n' );ylabel(yn);axis([0,n(end),min(xn),1.2*max(xn)]);实验二%时域采样理论验证程序exp2a.mTp=64/1000;%观察时间Tp=64微秒%产生M长采样序列x(n)% Fs=1000; T=1/Fs;Fs=1000; T=1/Fs;M=Tp*Fs; n=0:M-1;A=444.128; alph=pi*50*2^0.5; omega=pi*50*2^0.5;xnt=A*exp(-alph*n*T).*sin(omega*n*T);Xk=T*fft(xnt,M);%M点FFT[xnt)]yn='xa(nT)';subplot(3,2,1);tstem(xnt,yn); %调用自编绘图函数tstem绘制序列图box on; title('(a) Fs=1000Hz');k=0:M-1; fk=k/Tp;subplot(3,2,2); plot(fk,abs(Xk));title('(a) T*FT[xa(nT)],Fs=1000Hz');xlabel('f(Hz)'); ylabel('幅度');axis([0,Fs,0,1.2*max(abs(Xk))])%=================================% Fs=300Hz和Fs=200Hz的程序与上面Fs=1000Hz的程序完全相同。

数字信号处理实验全部程序MATLAB

实验一熟悉MATLAB环境一、实验目的(1)熟悉MATLAB的主要操作命令。

(2)学会简单的矩阵输入和数据读写。

(3)掌握简单的绘图命令。

(4)用MATLAB编程并学会创建函数。

(5)观察离散系统的频率响应。

二、实验内容认真阅读本章附录，在MATLAB环境下重新做一遍附录中的例子，体会各条命令的含义。

在熟悉了MATLAB基本命令的基础上，完成以下实验。

上机实验内容：(1)数组的加、减、乘、除和乘方运算。

输入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并用stem语句画出A、B、C、D、E、F、G。

实验程序：A=[1 2 3 4];B=[3 4 5 6];n=1:4;C=A+B;D=A-B;E=A.*B;F=A./B;G=A.^B;subplot(4,2,1);stem(n,A,'fill');xlabel ('时间序列n');ylabel('A');subplot(4,2,2);stem(n,B,'fill');xlabel ('时间序列n ');ylabel('B');subplot(4,2,3);stem(n,C,'fill');xlabel ('时间序列n ');ylabel('A+B');subplot(4,2,4);stem(n,D,'fill');xlabel ('时间序列n ');ylabel('A-B');subplot(4,2,5);stem(n,E,'fill');xlabel ('时间序列n ');ylabel('A.*B');subplot(4,2,6);stem(n,F,'fill');xlabel ('时间序列n ');ylabel('A./B');subplot(4,2,7);stem(n,G,'fill');xlabel ('时间序列n ');ylabel('A.^B');运行结果：(2)用MATLAB实现以下序列。

《数字信号处理》实验及代码

实验一用Matlab 进行数字信号处理一、实验目的1．掌握 MATLAB 基本的操作；2．学习典型的离散信号的Matlab 实现方法。

3．学习离散时间序列的基本运算：相加、相乘、移位等；二、实验内容1. 练习把y=sin(x), z=cos(x),u=2sin(x),v=sin(x)/cos(x)在[0，2pi]区间内的4个子图分别用不同的颜色、点型和线型绘制在同一个窗口中，并加上纵坐标、标题、图例和网格线。

2. 利用MATLAB 编程产生和绘制下列有限长序列： (1)单位脉冲序列()n δ (2)单位阶跃序列()u n (3)矩形序列8()R n(4)正弦序列()sin()53x n n ππ=+(5)实指数序列0.9()n u n3．上机调试并打印或记录实验结果。

4．完成实验报告。

三、实验结果1. 实验程序如下：x=0:pi/10:2*pi;y=sin(x); plot(x,y,'r*:');grid on %绘制网格线 hold on z=cos(x); plot(x,z,'y+-'); hold on u=2*sin(x); plot(x,u,'bx-'); hold on v=sin(x)/cos(x); plot(x,v,'ko-');hold onxlabel('x 轴'); %x轴注释ylabel('y 轴'); %y轴注释legend({'y=sin(x)','z=cos(x)', 'u=2sin(x) ', 'v=sin(x)/cos(x) '}); %图形注解2. 程序如下：n=-20:20;n0=0;n1=8;w0=pi/5; w1=pi/3;x1=[(n-n0)==0];x2=[(n-n0)>=0];x3=[(n-n0)>=0& (n-n1)<=0];x4=sin(w0*n+w1);x5=0.9.^n.*x2;subplot(511);stem(n,x1);axis([ -20 20 0 2]);ylabel('\sigma(n)'); subplot(512);stem(n,x2);axis([ -20 20 0 2]);ylabel('u(n)');subplot(513);stem(n,x3);axis([ -20 20 0 2]);ylabel('R8(n)'); subplot(514);stem(n,x4);axis([ -20 20 -2 2]);ylabel(' sin(w0n+w1) '); subplot(515);stem(n,x5);axis([ -20 20 0 2]);ylabel('0.9nu(n)'); xlabel('n');实验二离散信号与系统一、实验目的1．掌握卷积定理、熟悉离散信号和系统的时域特性；2．学习Matlab进行卷积计算方法；3．学习Matlab求解差分方程二、实验内容1．离散信号的基本运算：对序列x(n)={2,3,4,1,2,5} ,n=0,1,2,3,4,5,的移位、乘法、加法、翻转及尺度变换。

(完整版)数字信号处理试卷及答案_程培青(第三版),推荐文档

《数字信号处理》试卷 A 第 6 页（共 6 页）
数字信号处理基础试卷答案及评分标准
一、填空题：（共 28 分，每空 2 分）
7
建议收藏下载本文，以便随时学习！（1）2y(n)，y(n-3) （2）f≥2fs （3）N，抽样（4） X (k) xnWNnk n0
（5）递归型
（6）8
Z-1 0.5 -1.4
Z-1 -0.8 1
Z-1
Z-1
-0.8
1
3、
我去人也就有人！为UR扼腕入站内信不存在向你偶同意调剖沙龙课反倒是龙卷风前一天我分页符ZNBX吃噶十多个OK地方价格
复加所需时间T1 0.5106 N N 1 0.5106 512 511 0.130816s
所以T T1 T2 1.441536s
2、用 FFT 计算
复乘所需时间
T1
5 106
N 2
log2
N
5 106
512 2
log2
512
0.01152s
复加所需时间T2 0.5106 N log2 N 0.5106 512 log2 512 0.002304s
3、请画出 8 点的按频率抽取的（DIF）基-2 FFT 流图，要求输入自然数顺序，输出倒位序。
2、用级联型结构实现以下系统函数，试问一共能构成几种级联型网络，并画出结构图。
4Z 1Z 2 1.4Z 1 H (z) Z 0.5Z 2 0.9Z 0.8
专业班级：
学院名称
我去人也就有人！为UR扼腕入站内信不存在向你偶同意调剖沙龙课反倒是龙卷风前一天我分页符ZNBX吃噶十多个OK地方价格
。
A. 1
B.δ(w)
C. 2πδ(w)

数字信号处理上机实验及参考程序

数字信号处理上机实验及参考程序数字信号处理实验实验⼀离散时间信号与系统及MA TLAB实现1.单位冲激信号：n = -5:5;x = (n==0);subplot(122);stem(n, x);2.单位阶跃信号：x=zeros(1,11);n0=0;n1=-5;n2=5;n = n1:n2;x(:,n+6) = ((n-n0)>=0);stem(n,x);3.正弦序列：n = 0:1/3200:1/100;x=3*sin(200*pi*n+1.2);stem(n,x);4.指数序列n = 0:1/2:10;x1= 3*(0.7.^n);x2=3*exp((0.7+j*314)*n);subplot(221);stem(n,x1);subplot(222);stem(n,x2);5.信号延迟n=0:20;Y1=sin(100*n);Y2=sin(100*(n-3));subplot(221);stem(n,Y1);subplot(222);stem(n,Y2);6.信号相加X1=[2 0.5 0.9 1 0 0 0 0];X2=[0 0.1 0.2 0.3 0.4 0.5 0.6 0.7];X=X1+X2;stem(X);7.信号翻转X1=[2 0.5 0.9 1];n=1:4;X2=X1(5-n);subplot(221);stem(n,X1);subplot(222);stem(n,X2);8.⽤MATLAB计算序列{-2 0 1 –1 3}和序列{1 2 0 -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('幅度');9．⽤MA TLAB计算差分⽅程当输⼊序列为时的输出结果。

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('幅度')10.冲激响应impzN=64;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=impz(a,b,N);stem(k,y)xlabel('n');ylabel('幅度')11.传递函数频率响应a=[0.8 -0.44 0.36 0.22];%分⼦的系数数组b=[1 0.7 -0.45 -0.6];%分母的系数数组n=(0:500)*pi/500%在pi范围内取501个采样点[h,f]=freqz(a,b,n);%求系统的频率响应subplot(2,1,1),plot(n/pi,abs(h));grid%作系统的幅度频响图axis([0,1,1.1*min(abs(h)),1.1*max(abs(h))]);ylabel('幅度');subplot(2,1,2),plot(n/pi,angle(h));grid %作系统的相位频响图axis([0,1,1.1*min(angle(h)),1.1*max(angle(h))]);ylabel('相位');xlabel('以pi为单位的频率');12.系统零极点图a=[0.8 -0.44 0.36 0.22];b=[1 0.7 -0.45 -0.6];h=zplane(a,b);实验⼆离散信号变换1.解⽅程y(n)-2y(n-1)+3y(n-2)=4x(n)-5x(n-1)+6x(n-2)-7x(n-3)a=[4,-5,6,-7];b=[1,-2,3];n=[0:7]; x=ones(length(n));Y=[-1,1];X=[1,-1];xic=filtic(b,a,Y,X);y1=filter(b,a,x,xic)stem(n,y1);xlabel('n');ylabel('y(n)');2．对连续的单⼀频率周期信号按采样频率采样，截取长度N分别选N =20和N =16，观察其DFT结果的幅度谱。

数字信号处理实验1--5含代码

数字信号处理实验1--5含代码实验一离散时间信号的时域分析 1. 在MATLAB中利用逻辑关系式n,,0来实现序列，显示范围。

(产生如下，，,n,nn,n,n012图所示的单位脉冲信号的函数为impseq(n0,n1,n2)，程序如示例所示),3,n,10并利用impseq函数实现序列:; ，，，，，，yn,2,n,3，,n,6，，xn1nnnn120源代码:impseq.mfunction y=impseq(n0,n1,n2)n=[n1:n2]y=[(n-n0)==0]exp01-1.mfunction impseq(n0,n1,n2)n=-3:1:10y=2*impseq(3,-3,10)+impseq(6,-3,10);stem(n,y)n,,0，，2. 在MATLAB中利用逻辑关系式来实现序列，显示范围。

(自己编写un,nn,n,n012产生单位阶跃信号的函数，函数命名为stepseq(n0,n1,n2)) 并利用编写的stepseq函数实现序列: ，，，，，，yn,un，2，un,2,5,n,10源代码:stepseq.mfunction y=stepseq(n0,n1,n2)n=n1:1:n2y=[(n-n0)>=0]exp01-2.mfunction stepseq(n0,n1,n2)n=-5:1:20y=stepseq(-2,-5,20)+stepseq(2,-5,20)stem(n,y)3. 在MATLAB中利用数组运算符“.^”来实现一个实指数序列。

如: n ，，，，xn,0.30,n,15源代码:n=0:1:15;x=0.3.^nstem(n,x)4. 在MATLAB中调用函数sin或cos产生正余弦序列，如:π,, ，，，，xn,3sin0.4πn，，5cos0.3πn0,n,20,,5,,源代码:n=0:1:20x=11*sin(0.3*pi*n+pi/5)+5*cos(0.3*pi*n)stem(n,x)思考题:1.在MATLAB环境下产生单位脉冲序列和单位阶跃序列各有几种方法,如何使用,2.在MATLAB环境下进行序列的相乘运算时应注意什么问题,实验二离散时间系统的时域分析1. 在MATLAB中利用内部函数conv来计算两个有限长序列的卷积。

数字信号处理上机实验MATLAB程序及结果

研究生《数字信号处理》（dsp）上机MATLAB实验实验1熟悉MATLAB环境，会用FFT求信号频谱，相关方法去除噪声。

%生成一个带有噪声的正弦信号clear;t=linspace(-2*pi,2*pi,1000);y=sin(t)+rand(1,1000);plot(y);波形图：%把带有噪声的正弦信号进行FFT分解z=fft(y);plot(t,z);波形图：%用相关方法除噪声f=xcorr(y,y); plot(f);波形图：实验2利用相关算法求信号的幅值和相位%用相关性求信号的幅值与相位clear;t=linspace(-2*pi,2*pi,1000);s=5*sin(t+pi/6)+rand(1,1000);y=3*sin(t)+rand(1,1000);ss=s.*s;sy=s.*y;yy=y.*y;x1=length(ss);x2=length(yy);x3=length(sy);rx0=sum(ss)/x1;ry0=sum(yy)/x2;rxy0=sum(sy)/x3;A=sqrt(2*rx0);B=sqrt(2*ry0);C=acos(2*rxy0/(A*B));sprintf('A=%d,B=%d,C=%d',A,B,C)命令窗口输出结果：A=5.086632e+00,B=3.109529e+00,C=5.355640e-01 %用改进的相关算法求幅值和相位clear;t=linspace(-2*pi,2*pi,1000);s=5*sin(t+pi/6)+rand(1,1000);z=3*sin(t);z1=3*cos(t);zz=z.*z;sz=s.*z;sz1=s.*z1;x1=length(sz);x2=length(zz);x3=length(sz1);rz=sum(zz)/x2;rsz=sum(sz)/x1;rsz1=sum(sz1)/x3;B=sqrt(2*rz);C=atan(rsz1/rsz);A=2*rsz/(B*cos(C));sprintf('A=%d,B=%d,C=%d',A,B,C)命令窗口输出结果：A=4.986266e+00,B=2.998500e+00,C=5.264407e-01 实验三设计一个巴尔低通滤波器%设计一个巴尔沃特低通滤波器，fp=1khz,fs=2khz,ap=1db,as=40db;fp=1000;fs=2000;wp=2*pi*fp;ws=2*pi*fs;ap=1;as=40;[n,wc]=buttord(wp,ws,ap,as,'s');fprintf('阶数为%d\n',n);[num,den]=butter(n,wc,'s');%[z,p,k]=buttap(n);disp('分子多项式系数分别为');fprintf('%.4e\n',num);disp('分母多项式系数分别为');fprintf('%.4e\n',den);w=[wp ws];h1=freqs(num,den,w);fprintf('Ap=%.4e\n',-20*log10(abs(h1(1)))); fprintf('As=%.4e\n',-20*log10(abs(h1(2)))); w1=[0:200:12000*pi];h=freqs(num,den,w1);gain=20*log10(abs(h)); plot(w1/(2*pi),gain); xlabel('Hz');ylabel('dB');title('幅频特性');输出结果：阶数为8分子多项式系数分别为0.0000e+000.0000e+000.0000e+000.0000e+000.0000e+000.0000e+000.0000e+000.0000e+006.2187e+30分母多项式系数分别为1.0000e+003.6222e+046.5603e+087.7093e+126.4060e+163.8498e+201.6360e+244.5108e+276.2187e+30Ap=6.1668e-01As=4.0000e+01图形为实验四设计一个高通滤波器%设计一个高通滤波器fp=5khz,fs=1khz,ap<=1db,as>=40db; fp=5000;fs=1000;ap=1;as=40;wp=1/(2*pi*fp);ws=1/(2*pi*fs);[n,wc]=buttord(wp,ws,ap,as,'s');disp('阶数为');disp(n);[num,den]=butter(n,wc,'s');[numt,dent]=lp2hp(num,den,1);disp('高通滤波器分子多项式系数');fprintf('%.4e\n',numt)disp('高通滤波器分母多项式系数');fprintf('%.4e\n',dent)w=[1/wp 1/ws];h1=freqs(numt,dent,w);fprintf('Ap=%.4e\n',-20*log10(abs(h1(1))))fprintf('Bp=%.4e\n',-20*log10(abs(h1(2))))w1=[0:200:12000*pi];h=freqs(numt,dent,w1);gain=20*log10((abs(h)));plot(w1/(2*pi),gain);xlabel('频率');ylabel('幅值');输出结果：阶数为4高通滤波器分子多项式系数1.0000e+000.0000e+000.0000e+000.0000e+000.0000e+00高通滤波器分母多项式系数1.0000e+005.1920e+041.3478e+092.0497e+131.5585e+17Ap=1.0977e-01Bp=4.0000e+01图形实验5 用FFT设计低通滤波器%用fft变换设计低通滤波器clear;N=100;n=0:(N-1);fs=100;%采样频率f1=1;f2=10;y=sin(2*pi*f1/fs*n)+sin(2*pi*f2/fs*n);%原始信号subplot(211);plot(2*pi*f2/fs*n,y,'k');%显示波形xlabel('时间');ylabel('幅值');title('原始信号');ffty=fft(y);%进行快速傅里叶变换ffty1=ffty;ffty1(f2/(fs/N)+1)=0;%f2频率对应的幅值为零ffty1((N-f2)/(fs/N)+1)=0;%f2/2频率对应f2的对称频率值对应幅值为零yy=ifft(ffty1);%逆fft变换subplot(212);plot(2*pi*f2/fs*n,yy);xlabel('时间');title('滤波以后的低频信号');输出图形为实验六用FFT高通滤波器设计%用fft变换设计低通滤波器clear;N=100;n=0:(N-1);fs=100;%采样频率f1=1;f2=10;y=sin(2*pi*f1/fs*n)+sin(2*pi*f2/fs*n);%原始信号subplot(211);plot(2*pi*f2/fs*n,y,'k');%显示波形xlabel('时间');ylabel('幅值');title('原始信号');ffty=fft(y);%进行快速傅里叶变换ffty1=ffty;ffty1(f1/(fs/N)+1)=0;%f2频率对应的幅值为零ffty1((N-f1)/(fs/N)+1)=0;%f2/2频率对应f2的对称频率值对应幅值为零yy=ifft(ffty1);%逆fft变换subplot(212);plot(2*pi*f2/fs*n,yy);xlabel('时间');title('滤波以后的高频信号');输出图形为实验7 用双线性设计数字滤波器%用双线性设计数字滤波器,wp=0.2pi,ws=0.6Pi,ap=2db,as=15db; wp=0.2*pi;ws=0.6*pi;fs=1/2;ap=2;as=15;Wp=2*fs*tan(wp/2);Ws=2*fs*tan(ws/2);[n,wc]=buttord(Wp,Ws,ap,as,'s');[num,den]=butter(n,wc,'s');[numz,denz]=bilinear(num,den,fs);disp('分子多项式的系数');fprintf('%.4e\n',numz);disp('分母多项式的系数');fprintf('%.4e\n',denz);w=[wp ws];h=freqz(numz,denz,w);fprintf('ap=%.4f\n',-20*log10(abs(h(1))));fprintf('as=%.4f\n',-20*log10(abs(h(2)))); w=linspace(0,2*pi,1024);h=freqz(numz,denz,w);gain=20*log10(abs(h));plot(w/pi,gain);xlabel('Normalized frequency');ylabel('Gain,dB');axis([0 1 -50 0]);grid;输出结果分子多项式的系数1.5777e-013.1555e-011.5777e-01分母多项式的系数1.0000e+00-6.0620e-012.3730e-01ap=0.3945as=15.0000图形为实验8切比雪夫1型设计低通和高通滤波器%用切比雪夫1型设计低通滤波器%wp=2*pi*1000,ws=2*pi*2000,ap=1,as=40;clear;wp=2*pi*1000,ws=2*pi*2000,ap=1,as=40;[n,wc]=cheb1ord(wp,ws,ap,as,'s');[num,den]=cheby1(n,ap,wc,'s');disp('分子多项式的系数');fprintf('%.4f\n',num)disp('分母多项式的系数');fprintf('%.4f\n',den)w=[wp ws];h=freqs(num,den,w);fprintf('ap=%.4f\n',-20*log10(abs(h(1))));fprintf('as=%.4f\n',-20*log10(abs(h(2))));w=0:200:12000*pi;h=freqs(num,den,w);gain=20*log10(abs(h));plot(w/(2*pi),gain);输出结果分子多项式的系数0.00000.00000.00000.00000.00001202796127888840400.0000分母多项式的系数1.00005886.214566671782.4829241699153902.8650904788863911639.00001202796127888840400.0000ap=1.0000as=45.3060图形%用切比雪夫1型设计高通滤波器wp=2*pi*5000,ws=2*pi*1000,ap=1,as=40;wp=1/wp;ws=1/ws;[n,wc]=cheb1ord(wp,ws,ap,as,'s');[num,den]=cheby1(n,ap,wc,'s');[numz,denz]=lp2hp(num,den,1);w=[1/wp 1/ws];h=freqs(numz,denz,w);fprintf('ap=%.4e\n',-20*log10(abs(h(1)))); fprintf('as=%.4e\n',-20*log10(abs(h(2)))); w=0:200:12000*pi;h=freqs(numz,denz,w);gain=20*log10(abs(h));plot(w/(2*pi),gain);输出结果ap=1.0000e+00as=4.7847e+01输出图形实验9 用切比雪夫2型设计低通和高通滤波器%用切比雪夫2型设计低通滤波器%wp=2*pi*1000,ws=2*pi*2000,ap=1,as=40;clear;wp=2*pi*1000,ws=2*pi*2000,ap=1,as=40;[n,wc]=cheb2ord(wp,ws,ap,as,'s');[num,den]=cheby2(n,as,wc,'s');disp('分子多项式的系数');fprintf('%.4f\n',num)disp('分母多项式的系数');fprintf('%.4f\n',den)w=[wp ws];h=freqs(num,den,w);fprintf('ap=%.4f\n',-20*log10(abs(h(1))));fprintf('as=%.4f\n',-20*log10(abs(h(2))));w=0:200:12000*pi;h=freqs(num,den,w);gain=20*log10(abs(h)); plot(w/(2*pi),gain);输出结果分子多项式的系数0.0000566.39190.0000290688415499.92920.0000 29837909134910329000.0000分母多项式的系数1.000024344.5782296168842.7942 2252800523365.8032 10820843108280358.0000 29837909134910292000.0000 ap=1.0000as=44.1570图形%用切比雪夫1型设计高通滤波器wp=2*pi*5000;ws=2*pi*1000;ap=1;as=40;wp=1/wp;ws=1/ws;[n,wc]=cheb2ord(wp,ws,ap,as,'s');[num,den]=cheby2(n,as,wc,'s');[numz,denz]=lp2hp(num,den,1);w=[1/wp 1/ws];h=freqs(numz,denz,w);fprintf('ap=%.4e\n',-20*log10(abs(h(1))));fprintf('as=%.4e\n',-20*log10(abs(h(2))));w=0:200:12000*pi;h=freqs(numz,denz,w);gain=20*log10(abs(h));plot(w/(2*pi),gain);xlabel('频率');ylabel('幅值');title('幅频特性曲线');输出结果ap=1.0000e+00as=4.4779e+01输出图形实验10画一个信号的频谱图clear;fs=1000;f1=50;f2=100;f3=200;N=100t=0:N-1;data=25+22*sin(2*pi*f1/fs*t)+30*sin(2*pi*f2/fs*t)+50*sin(2*pi*f3/fs*t );%原始信号y = fft(data);yabs =abs(y);n = 0:N/2;yabs(1)=yabs(1)/2;stem(n*fs/N,2*yabs(1:N/2+1)/N);频谱图实验11 数字高通滤波器的设计%数字高通滤波器的设计，ws=0.2pi,wp=0.6pi,as=15,ap=2;用巴特沃斯滤波器设计。

数字信号处理(1) 代码

H=20*log10(abs(h)); %输出幅值转化为分贝单位
subplot(311); %把图像放到311位置
plot(t,x); %画出输入信号图像
[b,a]=butter(N,Wn,'s'); %巴特沃斯模拟滤波器设计
[b1,a1]=bilinear(b,a,0.5); %用双线性法，将模拟滤波器转换为数字滤波器
[h,omega]=freqz(b1,a1,256); %系统的绝对幅度响应（h包含了离散系统频响在 0~pi范围内256个频率等分点的值，omega则包含了范围内256个频率等分点）
ylabel('H (dB)'); %Y轴取值H
title('滤波器图像'); %图像题目
y=filter(b,a,x); %将信号进行滤波处理
title('滤波后的图像'); %图像题目
Rp=1; %幅度特性下降小于1DB
Rs=15; %幅度特性下降小于15DB
[N,Wn]=buttord(Wp,Ws,Rp,Rs,'s'); %巴特沃斯滤波器N估计
subplot(313); %把图像放到313位置
pl %网格布线
grid; %网格布线
title('输入信号'); %图像题目
subplot(312); %把图像放到312位置
x=sin(2*pi*20*t)+2*sin(2*pi*100*t)+5*sin(2*pi*200*t); %输入信号
Wp=tan(pi*30/Fs); %预畸变处理,Wp通带截至频率
Ws=tan(pi*60/Fs); %预畸变处理,Ws阻带截至频率

数字信号处理96142 - 道客巴巴

数字信号处理96142 - 道客巴巴信号处理一、实验要求编程利用FFT进行卷积计算，通过实验比较出快速卷积优越性。

写出实验报告。

二、实验原理利用FFT进行离散卷积的步骤归纳如下：⑴、设x(n)的列长为Np h(n)的列长为N2，要求 N-1N・1y(n)=x(n)*h(n)=自 x(k)h((n-k))NRN(n)=自 x伙)h(n0A-0 (2)、为使两有限长序列的线性卷积可用其圆周卷积來代替而不产生混淆，必须选择N2M+ N2-1o为使用基-2FFT来完成卷积计算，故要求22 (V是整数)。

用补零的办法使x(n), h(n)具有列长N,即 t x(n)n =0,12 技 x(n)= i fOn= Nl, Nl+1,技 NJ t h(n)n =0,1,2,技 N2八 h(n)= i, f fOn= N2, N2+1,技N-l ⑶为用圆周卷积定理计算线性卷积，先用FFT计算x(n) , h(n)的N点离散傅立叶变换 [2] x(n)揪揪瞅X(k)h(n)揪彌 H(k)(4)组成卷积 Y(k)=X(k)H(k)[4] (5)利用IFFT计算Y(k)的离散傅立叶逆变换得到线性卷积y(n)o 由于 N・IN・1 y(n)=a [(i/n)y伙)]wn虫咱[(1/川)厂(灯]\«令呵 *=0*=0 可见，y(n)可由求(l/N)Y*(k)的FFT再取共辄得到。

三、实验题目 (1)输入序列： dlllllllllllllmdllllllllllllmWIIIIIIIIIIII 严Flllllllllll““ 卷积输出: 正弦序列1 m (2)正弦序列与三角序列的卷积(正弦序列为两个周期，256点；三角序列为一个周期,256个点)输入序列： m (3)两个矩形序列的卷积(均为两个周期，256个点，占空比0、5)输入序列：卷积输出: 0 «：：»；：»>、；t、K、? •、 I » * A I 、i *150100500 矩形序列1150200 矩形序列21 020******* (4)单位冲击与正弦波(单位冲击序列为256个点，正弦序列为1、7个周期，256个点)输入序列：卷积输出: 正弦序列1 0、5 °?25030050100150200 卷积输出200・150・100-50050100150200250 任何序列与单位冲击序列的卷积为原序列，所以结果正确。

matlab上机程序(数字信号处理)

1x1=[1 1 1 1];x2=[1 2 3 4 4 3 2 1];n1=0:8;x3=cos(n1*pi/4);n2=0:8;x4=sin(n2*pi/8);figuresubplot(2,2,1) ;stem(x1);subplot(2,2,2) ;stem(x2);subplot(2,2,3) ;stem(x3);subplot(2,2,4) ;stem(x4);N1=8;F1x1=fft(x1,N1);F1x2=fft(x2,N1);F1x3=fft(x3,N1);F1x4=fft(x4,N1);subplot(2,2,1) ;stem(x1)grid on;title('juxing')xlabel('n')ylabel('x1(n)')figuresubplot(2,2,1) ;stem(abs(F1x1));subplot(2,2,2) ;stem(abs(F1x2));subplot(2,2,3) ;stem(abs(F1x3));subplot(2,2,4) ;stem(abs(F1x4));N3=256;F3x1=fft(x1,N3);F3x2=fft(x2,N3);F3x3=fft(x3,N3);F3x4=fft(x4,N3);w=2/N3*[0:N3-1];figuresubplot(2,2,1) ;plot(w,abs(F3x1));subplot(2,2,2) ;plot(w,abs(F3x2));subplot(2,2,3) ;plot(w,abs(F3x3));subplot(2,2,4) ;plot(w,abs(F3x4));N6=64;t=0:1/64:1/64*(N6-1);x6= cos(8*pi*t)+cos(16*pi*t)+cos(20*pi*t); figure;plot(x6)Fx6=fft(x6,N6);w6=2/N6*[0:N6-1];figure;plot(w6,abs(Fx6))grid on;title('juxing')xlabel('n')2--1clear all;close all;clcT=1;Wp=2/T*tan(0.2*pi/2);Ws=2/T*tan(0.3*pi/2);Rp=1;Rs=15;[N,Wc]=buttord(Wp,Ws,Rp,Rs,'s');[B,A]=butter(N,Wc,'s');[C,D]=bilinear(B,A,1/T);W=0:0.001*pi:0.5*pi;H=freqs(B,A,W);Hd=freqz(C,D,W);figuresubplot(3,1,1);plot(W/pi,abs(H))grid ontitle('模拟巴特沃斯滤波器')xlabel('Frequency/Hz')ylabel('Magnitude')subplot(3,1,2);plot(W/pi,abs(Hd))grid ontitle('数字巴特沃斯滤波器')xlabel('Didital Frequency/pi')ylabel('Magnitude')Hd_db=-20*log(abs(Hd(1)./(abs(Hd)+eps)));subplot(3,1,3);plot(W/pi,abs(Hd_db))grid ontitle('数字巴特沃斯滤波器波特图')xlabel('Didital Frequency/pi')ylabel('Magnitude')x=[-4,-2,0,-4,-6,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,-60,-84-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6, 4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0];M=512;fx=fft(x,M);f=1/M*[0:M-1]*250;y=filter(C,D,x);fy=fft(y,M)figuresubplot(2,2,1);plot(x);grid onsubplot(2,2,2);plot(f,abs(fx));grid on;subplot(2,2,3);plot(y);grid onsubplot(2,2,4);plot(f,abs(fy));grid on2----2clear all;close all;clcT=1;Wp=2/T*tan(0.2*pi/2);Ws=2/T*tan(0.3*pi/2);Rp=1;Rs=15;[N,Wc]=buttord(Wp,Ws,Rp,Rs,'s');[B,A]=butter(N,Wc,'s');[C,D]=bilinear(B,A,1/T);W=0:0.001*pi:0.5*pi;H=freqs(B,A,W);Hd=freqz(C,D,W);figuresubplot(3,1,1);plot(W/pi,abs(H))grid ontitle('模拟巴特沃斯滤波器')xlabel('Frequency/Hz')ylabel('Magnitude')figuresubplot(3,1,2);plot(W/pi,abs(Hd))grid ontitle('数字巴特沃斯滤波器')xlabel('Didital Frequency/pi')ylabel('Magnitude')Hd_db=-20*log(abs(Hd(1)./(abs(Hd)+eps)));figuresubplot(3,1,3);plot(W/pi,abs(Hd_db))grid ontitle('数字巴特沃斯滤波器波特图')xlabel('Didital Frequency/pi')ylabel('Magnitude')x=[-4,-2,0,-4,-6,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,-60,-84-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6, 4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0];M=512;fx=fft(x,M);y=filter(C,D,x);fy=fft(y,M);3function hd=ideal(N,wc)for n=0:N-1if n==(N-1)/2hd(n+1)=wc/pi;else hd(n+1)=sin(wc*(n-(N-1)/2))/(pi*(n-(N-1)/2));endendclear all;cla;close allN=33;wc=pi/4hd=ideal(N,wc);w1=boxcar(N);w2=hamming(N);w3=hann(N);w4=blackman(N);h1=hd.*w1';h2=hd.*w2';h3=hd.*w3';h4=hd.*w4';M=512;fh1=fft(h1,M);db1=-20*log10(abs(fh1(1)./(abs(fh1)+eps))); fh2=fft(h2,M);db2=-20*log10(abs(fh2(1)./(abs(fh2)+eps))); fh3=fft(h3,M);db3=-20*log10(abs(fh3(1)./(abs(fh3)+eps))); fh4=fft(h4,M);db4=-20*log10(abs(fh4(1)./(abs(fh4)+eps))); w=2/M*[0:M-1];figuresubplot(2,2,1);stem(h1)subplot(2,2,2);plot(w,abs(fh1))subplot(2,2,3);plot(w,db1)subplot(2,2,4);plot(w,angle(fh1))figuresubplot(2,2,1);stem(h2)subplot(2,2,2);plot(w,abs(fh2))subplot(2,2,3);plot(w,db2)subplot(2,2,4);plot(w,angle(fh2))figuresubplot(2,2,1);stem(h3)subplot(2,2,2);plot(w,abs(fh3))subplot(2,2,3);plot(w,db3)subplot(2,2,4);plot(w,angle(fh3))figuresubplot(2,2,1);stem(h4)subplot(2,2,2);plot(w,abs(fh4))subplot(2,2,3);plot(w,db4)subplot(2,2,4);plot(w,angle(fh4))。

数字信号处理实验(1-7)原始实验内容文档(含代码)

实验要求1.每个实验进行之前须充分预习准备，实验完成后一周内提交实验报告；2.填写实验报告时，分为实验题目、实验目的、实验内容、实验结果、实验小结五项；3.实验报告要求：实验题目、实验目的、实验内容、实验结果四项都可打印；但每次实验的实验内容中的重要代码（或关键函数）后面要用手工解释其作用。

实验小结必须手写！（针对以前同学书写实验报告时候抄写代码太费时间的现象，本期实验报告进行以上改革）。

实验一信号、系统及系统响应实验目的：1. 掌握使用MATLAB进行函数、子程序、文件编辑等基本操作；2. 编写一些数字信号处理中常用序列的3. 掌握函数调用的方法。

实验内容：1.在数字信号处理的基本理论和MATLAB信号处理工具箱函数的基础上，可以自己编写一些子程序以便调用。

（1）单位抽样序列δ(n-n0)的生成函数impseq.m（2）单位阶跃序列u(n-n0)的生成函数stepseq.m（3）两个信号相加的生成函数sigadd.m（4）两个信号相乘的生成函数sigmult.m（5）序列移位y(n)=x(n-n0)的生成函数sigshift.m（6）序列翻褶y(n)=x(-n)生成函数sigfold.m（7）奇偶综合函数evenodd.m（8）求卷积和2.产生系列序列，并绘出离散图。

（1） x1(n)=3δ(n-2)-δ(n+4) -5≤n≤5（2） x3(n)=cos(0.04πn)+0.2w(n) 0≤n≤50其中：w(n)是均值为0，方差为1 的白噪声序列。

3.设线性移不变系统的抽样响应h(n)=(0.9)^n*u(n),输入序列x(n)=u(n)-u(n-10),求系统的输出y(n).实验二系统响应及系统稳定性1.实验目的（1）掌握求系统响应的方法。

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

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

2.实验原理与方法在时域中，描写系统特性的方法是差分方程和单位脉冲响应，在频域可以用系统函数描述系统特性。

数字信号处理matlab pdf

数字信号处理matlab数字信号处理（Digital Signal Processing，DSP）是一种处理离散时间信号的方法，这些信号通常通过模拟信号进行数字化获得。

MATLAB是一个广泛使用的科学计算软件，可用于数字信号处理。

在MATLAB中进行数字信号处理的基本步骤通常包括：1.数据导入和预处理：MATLAB中可以方便地导入数字信号，包括音频、图像等。

然后可以对信号进行一些预处理操作，例如滤波、降噪等。

2.信号分析和特征提取：在信号预处理之后，可以进行更深入的分析，如频率分析、功率谱分析、相关性分析等。

此外，还可以提取信号的特征，例如频率、幅值、相位等。

3.信号处理算法实现：在MATLAB中，可以使用各种内置函数和工具箱来实现各种数字信号处理算法，如滤波器设计、频域变换、调制解调等。

4.结果可视化：MATLAB提供了强大的绘图和可视化工具，可以方便地显示信号处理的结果。

以下是一个简单的MATLAB代码示例，展示了如何读取一个音频文件并计算其功率谱：```matlab%读取音频文件[signal,fs]=audioread('filename.wav');%转换为单通道（如果需要）if size(signal,2)==2signal=sum(signal,2);end%计算功率谱[Pxx,F]=periodogram(signal,[],length(signal),fs);%绘制功率谱图figure;plot(F,10log10(Pxx/max(Pxx)));xlabel('Frequency(Hz)');ylabel('Power/Frequency(dB/Hz)');title('Power Spectrum');```请注意，这只是一个非常基础的示例。

实际应用中，数字信号处理可能涉及更复杂的算法和数据处理。

Matlab实现数字信号处理（源代码）

1 编制程序产生单位抽样序列δ(n) 及δ(n﹣学号)，并绘制出图形。

δ(n) %实现单位抽样序列δ(n)程序k=-30:30;delta=[zeros(1,30),1,zeros(1,30)];stem(k,delta)δ(n﹣学号)%实现单位抽样序列δ(n﹣22)程序k=-30:30;delta=[zeros(1,52),1,zeros(1,8)];stem(k,d2 编制程序产生单位阶跃序列u(n) 及u(n﹣学号)及u(n)﹣u(n﹣学号)，并绘制出图形。

%实现单位阶跃序列u(n)程序k=-30:30;delta=[zeros(1,30),1,ones(1,30)]; stem(k,delta)%实现阶跃序列u(n﹣22)程序k=-30:30;delta=[zeros(1,52),1,ones(1,8)];stem(k,delta)%实现阶跃序列u(n)﹣u(n﹣22)程序k=-30:30;(1,52),1,ones(1,8)];stem(k,delta)3 编制程序产生正弦序列x(n)= 学号sin(2πn)，x(n)= 学号sin(2πn/学号) 及x(n)= 学号sin(2n)，并绘制出图形。

%实现正弦序列x(n)= 22sin(2πn)程序n=0:100;x=22*sin(2*n*pi);stem(n,x)%实现正弦序列x(n)= 22sin(2πn/22)程序Array n=0:100;x=22*sin(2*n*pi/22);stem(n,x)n=0:100;x=22*sin(2*n);stem(n,x)4 编制程序产生复指数序列x(n)= e j学号n，并绘制图形。

%实现复指数序列x(n)= ej22n程序n=0:50;x=exp(i*22).^n;xr=real(x);xi=imag(x);xm=abs(x); xa=angle(x); figure;subplot(221);stem(n,xr);title('实部'); subplot(222);stem(n,xi);title('虚部'); subplot(223);stem(n,xm);title('模'); subplot(224);stem(n,xa);title('相角');实部204060虚部0204060模0204060相角5编制程序产生指数序列x(n)= a n，并绘制出图形。

数字信号处理课程设计(含完整代码)

课题一数字信号处理系统设计一、项目要求用本课程所学的数字信号处理理论知识，设计一个具有信号的采集、处理、传输、显示和存储等功能的系统，内容如下：１、录制一段语音信号，并对录制的语音信号进行采样（采样频率可取fs=22050Hz）;２、画出采样后的语音信号的时域波形和频谱图；３、滤波器的性能指标：低通滤波器：通带边界频率ｆｐ＝１ｋＨｚ，通带最大衰减Ａｐ＝１ｄB；阻带边界频率ｆｐ＝１.2ｋＨｚ，阻带最小衰减Ａｐ＝１00ｄB；高通滤波器：通带边界频率ｆｐ＝5ｋＨｚ，通带最大衰减Ａｐ＝１ｄB；阻带边界频率ｆｐ＝4.8ｋＨｚ，阻带最小衰减Ａｐ＝１00ｄB；带通滤波器：通带上限截止频率fp2=3kHz, 通带下限截止频率fp1=1.2kHz;阻带上限截止频率fs2=3.2kHz, 通带下限截止频率fp1=1kHz;通带最大衰减Ap=1dB, 阻带最小衰减As=100dB;采用双线性变换法设计滤波器，并画出滤波器的频率响应；４、用自己设计的滤波器对采样的信号进行滤波，画出滤波后信号的时域波形和频谱图，并对滤波前后的信号进行对比，分析信号的变化；５、回放语音６、用ＧＵＩ设计一个信号系统的用户界面。

二、实验所要用到的ＭＡＴＬＡＢ函数１、语音信号的采样与播放wavread();[y，fs，bite]=wavread（）；%语音信号的采样sound（y，fs，bite）；%播放语音2.滤波器：IIR：butte（）；%巴特沃思滤波器cheby1(); %切比雪夫I滤波器elliptical（）；%椭圆滤波器3.频率响应：[h，f]=freqz（b,a,n,fs）freqz(b,a,n,fs)5.快速傅里叶变换fft (x, n)6.画曲线plot(x, y)stem(x, y)7.在MATLAB中，设计辅助低通原型巴特沃思和切比雪夫滤波器的阶数和截止频率；1）利用buttord和cheblord确定阶数；2）[num，den]=butter（N，Wn），[num，den]=cheby1（N，Wn）3）lp2hp,lp2bp,lp2bs可以完成低通滤波器到高通，带通，带阻的转换4）使用biliner对模拟滤波器进行双线性变换，求得数字滤波器的传输函数系数三、数字滤波器（代码）1.IIR低通滤波器fp=1000;ft=5500;fs=1200;wp=2*pi*fp/ft;ws=2*pi*fs/ft;Fp=2*ft*tan(wp/2);Fs=2*ft*tan(ws/2);[n,Omgc]=buttord(Fp,Fs,1,100,'s');[z,p,k]=buttap(n);B=k*real(poly(z));A=real(poly(p));[b1,a1]=lp2lp(B,A,Omgc);[ba1,aa1]=bilinear(b1,a1,ft); %模拟转数字[Ha,w]=freqz(ba1,aa1); %求频率相应plot(w*ft/(2*pi),20*log10(abs(Ha)));title('IIR 低通滤波器');xlabel('频率/HZ');ylabel('幅值');2.IIR高通滤波器fp=5000;ft=25000;%取抽样频率•fs=4800;Rp=1;As=100;wp=2*pi*fp/ft;%通带频率ws=2*pi*fs/ft;%截止频率Fp=2*ft*tan(wp/2);Fs=2*ft*tan(ws/2);[n,Omgc]=ellipord(Fp,Fs,Rp,As,'s'); %计算阶数n和截止频率[z,p,k]=ellipap(n,Rp,As);B3=k*real(poly(z));A3=real(poly(p));[b3,a3]=lp2hp(B3,A3,Omgc);[ba3,aa3]=bilinear(b3,a3,ft);[Ha,w]=freqz(ba3,aa3);plot(w*ft/(2*pi),20*log10(abs(Ha)));title('IIR高通滤波器');xlabel('频率/HZ');ylabel('幅值');3.IIR带通滤波器fp1=1200;fp2=3000;fs1=1000;fs2=3200;ft=10000; As=100;Rp=1;wp1=2*pi*fp1/ft;wp2=2*pi*fp2/ft;ws1=2*pi*fs1/ft;ws2=2*pi*fs2/ft;Fp1=2*ft*tan(wp1/2);Fp2=2*ft*tan(wp2/2);Fp=[Fp1,Fp2];Fs1=2*ft*tan(ws1/2);Fs2=2*ft*tan(ws2/2);Fs=[Fs1,Fs2];bw=Fp2-Fp1;w0=sqrt(Fp1*Fp2);%通带宽和中心频率[n,Omgn]=cheb1ord(Fp,Fs,Rp,As,'s');[z,p,k]=cheb1ap(n,Rp);B2=k*real(poly(z));A2=real(poly(p));[b2,a2]=lp2bp(B2,A2,w0,bw);[ba2,aa2]=bilinear(b2,a2,ft);[Ha,w]=freqz(ba2,aa2);plot(w*ft/(2*pi),20*log10(abs(Ha)));title('IIR 带通滤波器'); xlabel('频率/HZ'); ylabel('幅值');四、总代码function varargout = bzh(varargin)% BZH M-file for bzh.fig% BZH, by itself, creates a new BZH or raises the existing% singleton*.%% H = BZH returns the handle to a new BZH or the handle to% the existing singleton*.%% BZH('CALLBACK',hObject,eventData,handles,...) calls the local% function named CALLBACK in BZH.M with the given input arguments.%% BZH('Property','Value',...) creates a new BZH or raises the % existing singleton*. Starting from the left, property value pairs are% applied to the GUI before ko_OpeningFunction gets called. An % unrecognized property name or invalid value makes property application% stop. All inputs are passed to bzh_OpeningFcn via varargin. %% *See GUI Options on GUIDE's Tools menu. Choose "GUI allowsonly one% instance to run (singleton)".%% See also: GUIDE, GUIDATA, GUIHANDLES% Copyright 2002-2003 The MathWorks, Inc.% Edit the above text to modify the response to help bzh% Last Modified by GUIDE v2.5 09-Jan-2014 08:54:22% Begin initialization code - DO NOT EDITgui_Singleton = 1;gui_State = struct('gui_Name', mfilename, ...'gui_Singleton', gui_Singleton, ...'gui_OpeningFcn', @bzh_OpeningFcn, ...'gui_OutputFcn', @bzh_OutputFcn, ...'gui_LayoutFcn', [] , ...'gui_Callback', []);if nargin && ischar(varargin{1})gui_State.gui_Callback = str2func(varargin{1});endif nargout[varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:}); elsegui_mainfcn(gui_State, varargin{:});end% End initialization code - DO NOT EDIT% --- Executes just before bzh is made visible.function bzh_OpeningFcn(hObject, eventdata, handles, varargin)% This function has no output args, see OutputFcn.% hObject handle to figure% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA)% varargin command line arguments to bzh (see VARARGIN)% Choose default command line output for bzhhandles.output = hObject;% Update handles structureguidata(hObject, handles);% UIWAIT makes bzh wait for user response (see UIRESUME)% uiwait(handles.figure1);% --- Outputs from this function are returned to the command line. function varargout = bzh_OutputFcn(hObject, eventdata, handles) % varargout cell array for returning output args (see VARARGOUT); % hObject handle to figure% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA)% Get default command line output from handles structurevarargout{1} = handles.output;% --- Executes on button press in tag_start.function tag_start_Callback(hObject, eventdata, handles)% hObject handle to tag_start (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) axes(handles.axes1)cla;global z0;global yy;fs=22050;nbits=32;[yy,fs,nbits]=wavread('D:\f.wav'); %语音信号加载sound(yy,fs); %回放语音z0=yy;axes(handles.axes1);plot(yy);title('时域采样信号波形');grid onn=length(yy); %求出语音信号的长度Y=fft(yy,n); %傅里叶变换axes(handles.axes2);plot(20*log10(abs(Y)));title('时域采样信号频谱');guidata(hObject,handles);grid on% --- Executes on button press in tag_ditong.function tag_ditong_Callback(hObject, eventdata, handles)% hObject handle to tag_ditong (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) axes(handles.axes1)cla;global z1;global yy;global ba1;global aa1;global Fp;global Fs;fp=1000;ft=5500;fs=1200;wp=2*pi*fp/ft;ws=2*pi*fs/ft;Fp=2*ft*tan(wp/2);Fs=2*ft*tan(ws/2);[n,Omgc]=buttord(Fp,Fs,1,100,'s');[z,p,k]=buttap(n);B=k*real(poly(z));A=real(poly(p));[b1,a1]=lp2lp(B,A,Omgc);[ba1,aa1]=bilinear(b1,a1,ft); %模拟转数字z1=filter(ba1,aa1,yy);[Ha,w]=freqz(ba1,aa1); %求频率相应plot(w*ft/(2*pi),20*log10(abs(Ha)));title('IIR 低通滤波器');xlabel('频率/HZ');ylabel('幅值');% --- Executes on button press in tag_gaotong.function tag_gaotong_Callback(hObject, eventdata, handles)% hObject handle to tag_gaotong (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) axes(handles.axes2)cla;global yy;global z2;global B3;global ba3;global aa3;global Fp;global Fs;fp=5000;ft=25000;%取抽样频率•fs=4800;Rp=1;As=100;wp=2*pi*fp/ft;%通带频率ws=2*pi*fs/ft;%截止频率Fp=2*ft*tan(wp/2);Fs=2*ft*tan(ws/2);[n,Omgc]=ellipord(Fp,Fs,Rp,As,'s'); %计算阶数n和截止频率[z,p,k]=ellipap(n,Rp,As);B3=k*real(poly(z));A3=real(poly(p));[b3,a3]=lp2hp(B3,A3,Omgc);[ba3,aa3]=bilinear(b3,a3,ft);z2=filter(ba3,aa3,yy);[Ha,w]=freqz(ba3,aa3);plot(w*ft/(2*pi),20*log10(abs(Ha)));title('IIR高通滤波器');xlabel('频率/HZ');ylabel('幅值');% --- Executes on button press in tag_daitong.function tag_daitong_Callback(hObject, eventdata, handles)% hObject handle to tag_daitong (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) axes(handles.axes3)cla;global yy;global z3;global B2;global ba2;global aa2;fp1=1200;fp2=3000;fs1=1000;fs2=3200;ft=10000;As=100;Rp=1;wp1=2*pi*fp1/ft;wp2=2*pi*fp2/ft;ws1=2*pi*fs1/ft;ws2=2*pi*fs2/ft;Fp1=2*ft*tan(wp1/2);Fp2=2*ft*tan(wp2/2);Fp=[Fp1,Fp2];Fs1=2*ft*tan(ws1/2);Fs2=2*ft*tan(ws2/2);Fs=[Fs1,Fs2];bw=Fp2-Fp1;w0=sqrt(Fp1*Fp2);%通带宽和中心频率[n,Omgn]=cheb1ord(Fp,Fs,Rp,As,'s');[z,p,k]=cheb1ap(n,Rp);B2=k*real(poly(z));A2=real(poly(p));[b2,a2]=lp2bp(B2,A2,w0,bw);[ba2,aa2]=bilinear(b2,a2,ft);z3=filter(ba2,aa2,yy);[Ha,w]=freqz(ba2,aa2);plot(w*ft/(2*pi),20*log10(abs(Ha)));title('IIR 带通滤波器');xlabel('频率/HZ');ylabel('幅值');% --- Executes on button press in tag_dibo.function tag_dibo_Callback(hObject, eventdata, handles)% hObject handle to tag_dibo (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) axes(handles.axes3)cla;global yy;global z1;global ba1;global aa1;m4=fft(z1);plot(z1,'y');title('IIR滤波后的信号波形');xlabel('时间/t');ylabel('幅值');axes(handles.axes4)cla;plot(abs(m4),'y');title('IIR滤波后信号频谱');xlabel('频率/HZ');ylabel('幅值');% --- Executes on button press in tag_gaobo.function tag_gaobo_Callback(hObject, eventdata, handles)% hObject handle to tag_gaobo (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) axes(handles.axes3)cla;global yy;global z2;global B3;global ba3;global aa3;m5=fft(z2);plot(z2,'y');title('IIR滤波后的信号波形');xlabel('时间/t');ylabel('幅值');axes(handles.axes4)cla;plot(abs(m5),'y');title('IIR滤波后信号频谱');xlabel('频率/HZ');ylabel('幅值');% --- Executes on button press in tag_daibo.function tag_daibo_Callback(hObject, eventdata, handles)% hObject handle to tag_daibo (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) axes(handles.axes3)cla;global yy;global z3;global B2;global ba2;global aa2;m6=fft(z3);plot(z3,'y');title('IIR滤波后的信号波形');xlabel('时间/t');ylabel('幅值');axes(handles.axes4)cla;plot(abs(m6),'y');title('IIR滤波后信号频谱');xlabel('频率/HZ');ylabel('幅值');% --- Executes on selection change in tag_choose.function tag_choose_Callback(hObject, eventdata, handles)% hObject handle to tag_choose (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA)% Hints: contents = get(hObject,'String') returns tag_choose contents as cell array% contents{get(hObject,'Value')} returns selected item from tag_choose% --- Executes during object creation, after setting all properties. function tag_choose_CreateFcn(hObject, eventdata, handles)% hObject handle to tag_choose (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB % handles empty - handles not created until after all CreateFcns called% Hint: popupmenu controls usually have a white background on Windows. % See ISPC and COMPUTER.if ispcset(hObject,'BackgroundColor','white');elseset(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundCo lor'));end% --- Executes on button press in tag_return.function tag_return_Callback(hObject, eventdata, handles)% hObject handle to tag_return (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) global z0;global z1;global z2;global z3;global yy;val=get(handles.tag_choose,'Value');if(val==1)[yy,fs,nbits]=wavread('D:\f.wav'); %语音信号加载sound(yy,fs); %回放语音endif(val==2)sound(z1,22050);endif(val==3)sound(z3,22050);endif(val==4)sound(z2,22050);end% --- Executes on mouse press over axes background.function axes1_ButtonDownFcn(hObject, eventdata, handles)% hObject handle to axes1 (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA)% --- If Enable == 'on', executes on mouse press in 5 pixel border. % --- Otherwise, executes on mouse press in 5 pixel border or over tag_start.function tag_start_ButtonDownFcn(hObject, eventdata, handles)% hObject handle to tag_start (see GCBO)% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA)采样后信号滤波器数字波形IIR低通滤波后波形图IIR高通滤波后波形图IIR带通滤波后波形图。

数字信号处理编程一

编程作业实验：自己构造含高频噪声和低频有用信号的复杂信号，分析信号的频谱特征，据此给定一个低通滤波器的指标，分别借助于Butterworth、切比雪夫1型、切比雪夫2型、椭圆模拟低通滤波器，用MATLAB设计此低通数字滤波器并滤波。

实验要求：1)给出完整的程序；2)画出幅度谱和对应的增益图（dB）；3)分析不同方式所需的滤波器阶数；4)画出并分析每种滤波器的滤波前后的结果；实验过程：1高频噪声信号和低频有用信号程序：x=0:100;y=sin(0.1*pi*x);noise=sin(0.5*pi*x);y1=y+noise;figureplot(x,y,'r',x,y1,'b');title('红色的为原始信号，蓝色为加噪声后的信号');figurestem(x/100,abs(fft(y1)));xlabel('\omega/\pi');ylabel('幅度');title('加了噪声信号的频谱');运行结果图：0102030405060708090100-2-1.5-1-0.50.511.52红色的为原始信号，蓝色为加噪声后的信号ω/π幅度加了噪声信号的频谱结果分析：这次有用信号是y=sin(0.1*pi*x);频率为0.1π；噪声信号是noise=sin(0.5*pi*x);频率是0.5π；因此设置边界频率为wp=0.1π，通带的纹波不超过0.5dB，在阻带边界的频率为ws=0.5π处的最小阻带衰减为15dB；2、Butterworth滤波器的设计程序：x=0:100;y=sin(0.1*pi*x);noise=sin(0.5*pi*x);y1=y+noise;%自己设置的巴特沃斯滤波器的指标wp=0.1;ws=0.5;rp=0.5;rs=15;%构造的巴特沃斯滤波器[n,wn]=buttord(wp,ws,rp,rs);[a,b]=butter(n,wn);[h,omega] = freqz(a,b,256);gain =20*log10(abs(h));figureplot (omega/pi,gain);grid;xlabel('\omega/\pi'); ylabel('幅度, dB');title('构造的巴特沃斯滤波器');lvbo=filtfilt(a,b,y1);figureplot(x,y,'r',x,lvbo,'b');title('红色的为原始信号，蓝色为用巴特沃斯滤波器过滤的信号'); disp('巴特沃斯滤波器所需要的阶数n=');n程序运行结果：00.10.20.30.40.50.60.70.80.91-120-100-80-60-40-20ω/π幅度, d B 构造的巴特沃斯滤波器实验结果分析：巴特沃斯滤波器所需的阶数是n=33、切比雪夫1型滤波器程序：x=0:100;y=sin(0.1*pi*x);noise=sin(0.5*pi*x);y1=y+noise;%自己设置的切比雪夫1型波器的指标wp=0.1;ws=0.5;0102030405060708090100-1.5-1-0.50.511.5红色的为原始信号，蓝色为用巴特沃斯滤波器过滤的信号rp=0.5;rs=15;%构造的切比雪夫1型滤波器[n,wn]=cheb1ord(wp,ws,rp,rs);[a,b] = cheby1(n,rp,wn);[h,omega] = freqz(a,b,256);gain =20*log10(abs(h));figureplot (omega/pi,gain);grid;xlabel('\omega/\pi'); ylabel('幅度, dB');title('构造的切比雪夫1型滤波器');lvbo=filtfilt(a,b,y1);figureplot(x,y,'r',x,lvbo,'b');title('红色的为原始信号，蓝色为用切比雪夫1型波器过滤的信号');disp('切比雪夫1型所需要的阶数n=');nn运行结果图：00.10.20.30.40.50.60.70.80.91-120-100-80-60-40-20ω/π幅度, d B 构造的切比雪夫1型滤波器0102030405060708090100-1-0.8-0.6-0.4-0.20.20.40.60.81红色的为原始信号，蓝色为用切比雪夫1型波器过滤的信号运行结果分析：切比雪夫1型滤波器所需的阶数是2阶3、切比雪夫2型滤波器程序：x=0:100;y=sin(0.1*pi*x);noise=sin(0.5*pi*x);y1=y+noise;%自己设置的切比雪夫2型波器的指标wp=0.1;ws=0.5;rp=0.5;rs=15;%构造的切比雪夫2型滤波器[n,wn]=cheb2ord(wp,ws,rp,rs);[a,b] = cheby2(n,rs,wn);[h,omega] = freqz(a,b,256);gain =20*log10(abs(h));figureplot (omega/pi,gain);grid;xlabel('\omega/\pi'); ylabel('幅度, dB');title('构造的切比雪夫2型滤波器');lvbo=filtfilt(a,b,y1);figureplot(x,y,'r',x,lvbo,'b');title('红色的为原始信号，蓝色为用切比雪夫2型波器过滤的信号');disp('切比雪夫2型所需要的阶数n=');nn运行结果图：00.10.20.30.40.50.60.70.80.91-60-50-40-30-20-10ω/π幅度, d B 构造的切比雪夫2型滤波器0102030405060708090100-1.5-1-0.50.511.5红色的为原始信号，蓝色为用切比雪夫2型波器过滤的信号运行结果分析：切比雪夫2型滤波器所需的阶数是2阶4、椭圆滤波器程序：x=0:200;y=sin(0.1*pi*x);noise=sin(0.5*pi*x);y1=y+noise;%自己设置的椭圆波器的指标wp=0.1;ws=0.5;rp=0.5;rs=15;%构造的椭圆滤波器[n,wn]=ellipord(wp,ws,rp,rs);[a,b] = ellip(n,rp,rs,wn);[h,omega] = freqz(a,b,256);gain =20*log10(abs(h));figureplot (omega/pi,gain);grid;xlabel('\omega/\pi'); ylabel('幅度, dB');title('构造的椭圆滤波器');lvbo=filtfilt(a,b,y1);figureplot(x,y,'r',x,lvbo,'b');title('红色的为原始信号，蓝色为用切比雪夫2型波器过滤的信号'); disp('椭圆所需要的阶数n=');n程序运行图：00.10.20.30.40.50.60.70.80.91-60-50-40-30-20-10ω/π幅度, d B 构造的椭圆滤波器运行结果分析：椭圆滤波器所需要的阶数是2阶实验结果分析：各个滤波器滤波所需的阶数四种滤波器没有显示出太大的差别。

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

%%%%%%%%%%%% fs=300HzTp=64/1000;%观察时间Tp=64微秒%产生M长采样序列x(n)% Fs=1000; T=1/Fs;Fs=300; T=1/Fs;M=Tp*Fs; n=0:M-1;A=444.128; alph=pi*50*2^0.5; omega=pi*50*2^0.5;xnt=A*exp(-alph*n*T).*sin(omega*n*T);Xk=T*fft(xnt,M);%M点FFT[xnt)]yn='xa(nT)';figure(2);subplot(3,2,1);tstem(xnt,yn); %调用自编绘图函数tstem绘制序列图box on; title('(b) Fs=300Hz');k=0:M-1; fk=k/Tp;subplot(3,2,2); plot(fk,abs(Xk));title('(b) T*FT[xa(nT)],Fs=1000Hz');xlabel('f(Hz)'); ylabel('幅度');axis([0,Fs,0,1.2*max(abs(Xk))])%%%%%%%%%%%%%%%%%%% fs=200HzTp=64/1000;%观察时间Tp=64微秒%产生M长采样序列x(n)% Fs=1000; T=1/Fs;Fs=200; T=1/Fs;M=Tp*Fs; n=0:M-1;A=444.128; alph=pi*50*2^0.5; omega=pi*50*2^0.5;xnt=A*exp(-alph*n*T).*sin(omega*n*T);Xk=T*fft(xnt,M);%M点FFT[xnt)]yn='xa(nT)';figure(3);subplot(3,2,1);tstem(xnt,yn); %调用自编绘图函数tstem绘制序列图box on; title('(c) Fs=200Hz');k=0:M-1; fk=k/Tp;subplot(3,2,2); plot(fk,abs(Xk));title('(c) T*FT[xa(nT)],Fs=1000Hz');xlabel('f(Hz)'); ylabel('幅度');axis([0,Fs,0,1.2*max(abs(Xk))])%频域采样理论验证程序exp2b.mM=27;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,1024);%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:1023;wk=2*k/1024;%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])function tstem(xn,yn)n = 0:length(xn)-1;stem(n,xn,'.');xlabel('n');ylabel('yn');%xlabel('n' );ylabel(yn);axis([0,n(end),min(xn),1.2*max(xn)]);实验三%实验三程序exp3.m%用FFT对信号作频谱分析%clear all;close all%实验内容(1)================================ x1n=[ones(1,4)];%产生序列向量x1(n)=R4(n)M=8;xa=1:(M/2); xb=(M/2):-1:1;x2n=[xa,xb]; %产生长度为8的三角波序列x2(n) x3n=[xb,xa];X1k8=fft(x1n,8); %计算x1n的8点DFTX1k16=fft(x1n,16); %计算x1n的16点DFTX2k8=fft(x2n,8); %计算x1n的8点DFTX2k16=fft(x2n,16); %计算x1n的16点DFTX3k8=fft(x3n,8); %计算x1n的8点DFTX3k16=fft(x3n,16); %计算x1n的16点DFT%以下绘制幅频特性曲线subplot(2,2,1);mstem(X1k8);%绘制8点DFT的幅频特性图title('(1a) 8点DFT[x_1(n)]');xlabel('ω/π');ylabel('幅度');axis([0,2,0,1.2*max(abs(X1k8))])subplot(2,2,3);mstem(X1k16);%绘制16点DFT的幅频特性图title('(1b)16点DFT[x_1(n)]');xlabel('ω/π');ylabel('幅度');axis([0,2,0,1.2*max(abs(X1k16))])figure(2)subplot(2,2,1);mstem(X2k8);%绘制8点DFT的幅频特性图title('(2a) 8点DFT[x_2(n)]');xlabel('ω/π');ylabel('幅度');axis([0,2,0,1.2*max(abs(X2k8))])subplot(2,2,2);mstem(X2k16);%绘制16点DFT的幅频特性图title('(2b)16点DFT[x_2(n)]');xlabel('ω/π');ylabel('幅度');axis([0,2,0,1.2*max(abs(X2k16))])subplot(2,2,3);mstem(X3k8);%绘制8点DFT的幅频特性图title('(3a) 8点DFT[x_3(n)]');xlabel('ω/π');ylabel('幅度');axis([0,2,0,1.2*max(abs(X3k8))])subplot(2,2,4);mstem(X3k16);%绘制16点DFT的幅频特性图title('(3b)16点DFT[x_3(n)]');xlabel('ω/π');ylabel('幅度');axis([0,2,0,1.2*max(abs(X3k16))])%实验内容(2) 周期序列谱分析=================== N=8;n=0:N-1; %FFT的变换区间N=8x4n=cos(pi*n/4);x5n=cos(pi*n/4)+cos(pi*n/8);X4k8=fft(x4n); %计算x4n的8点DFTX5k8=fft(x5n); %计算x5n的8点DFTN=16;n=0:N-1; %FFT的变换区间N=16x4n=cos(pi*n/4);x5n=cos(pi*n/4)+cos(pi*n/8);X4k16=fft(x4n); %计算x4n的16点DFTX5k16=fft(x5n); %计算x5n的16点DFTfigure(3)subplot(2,2,1);mstem(X4k8);%绘制8点DFT的幅频特性图title('(4a) 8点DFT[x_4(n)]');xlabel('ω/π');ylabel('幅度');axis([0,2,0,1.2*max(abs(X4k8))])subplot(2,2,3);mstem(X4k16);%绘制16点DFT的幅频特性图title('(4b)16点DFT[x_4(n)]');xlabel('ω/π');ylabel('幅度');axis([0,2,0,1.2*max(abs(X4k16))])subplot(2,2,2);mstem(X5k8);%绘制8点DFT的幅频特性图title('(5a) 8点DFT[x_5(n)]');xlabel('ω/π');ylabel('幅度');axis([0,2,0,1.2*max(abs(X5k8))])subplot(2,2,4);mstem(X5k16);%绘制16点DFT的幅频特性图title('(5b)16点DFT[x_5(n)]');xlabel('ω/π');ylabel('幅度');axis([0,2,0,1.2*max(abs(X5k16))])%实验内容(3) 模拟周期信号谱分析================= figure(4)Fs=64;T=1/Fs;N=16;n=0:N-1; %FFT的变换区间N=16x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T);%对x6(t)16点采样X6k16=fft(x6nT); %计算x6nT的16点DFTX6k16=fftshift(X6k16); %将零频率移到频谱中心Tp=N*T;F=1/Tp; %频率分辨率Fk=-N/2:N/2-1;fk=k*F;%产生16点DFT对应的采样点频率（以零频率为中心）subplot(3,1,1);stem(fk,abs(X6k16),'.');box on %绘制8点DFT的幅频特性图title('(6a) 16点|DFT[x_6(nT)]|');xlabel('f(Hz)'); %%%%%%%%%%%%%%%%%%ylabel('幅度');axis([-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k16))])N=32;n=0:N-1; %FFT的变换区间N=16x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T);%对x6(t)32点采样X6k32=fft(x6nT); %计算x6nT的32点DFTX6k32=fftshift(X6k32); %将零频率移到频谱中心Tp=N*T;F=1/Tp; %频率分辨率Fk=-N/2:N/2-1;fk=k*F;%产生16点DFT对应的采样点频率（以零频率为中心）subplot(3, 1, 2);stem(fk, abs(X6k32), '.');box on %绘制8点DFT的幅频特性图title('(6b) 32点|DFT[x_6(nT)]|');xlabel('f(Hz)');ylabel('幅度');axis([-N*F/2-1, N*F/2-1,0, 1.2*max(abs(X6k32))])N=64;n=0:N-1; %FFT的变换区间N=16x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T);%对x6(t)64点采样X6k64=fft(x6nT); %计算x6nT的64点DFTX6k64=fftshift(X6k64); %将零频率移到频谱中心Tp=N*T;F=1/Tp; %频率分辨率Fk=-N/2:N/2-1;fk=k*F;%产生16点DFT对应的采样点频率（以零频率为中心）subplot(3,1,3);stem(fk,abs(X6k64),'.');box on%绘制8点DFT的幅频特性图title('(6a) 64点|DFT[x_6(nT)]|');xlabel('f(Hz)');ylabel('幅度');axis([-N*F/2-1,N*F/2-1,0, 1.2*max(abs(X6k64))])function mstem(Xk)%mstem(Xk)绘制频域采样序列向量Xk的幅频特性图M=length(Xk);k=0:M-1;wk=2*k/M;%产生M点DFT对应的采样点频率（关于pi归一化值) stem(wk,abs(Xk),'.');box on;%绘制M点DFT的幅频特性图?xlabel('w/\pi');ylabel('幅度');axis([0,2,0,1.2*max(abs(Xk))]);实验四%实验四程序exp4.m% IIR数字滤波器设计及软件实现%clear all; close allFs=10000; T=1/Fs; %采样频率%调用信号产生函数mstg产生由三路抑制载波调幅信号相加构成的复合信号st st=mstg;%低通滤波器设计与实现=========================fp=280; fs=450;wp=2*fp/Fs; ws=2*fs/Fs; rp=0.1; rs=60;%DF指标（低通滤波器的通、阻带边界频率）[N,wp]=ellipord(wp,ws,rp,rs);%调用ellipord计算椭圆DF阶数N和通带截止频率wp[B,A]=ellip(N,rp,rs,wp);%调用ellip计算椭圆带通DF系统函数系数向量B和Ay1t=filter(B,A,st); %滤波器软件实现% 低通滤波器设计与实现绘图部分figure(2); subplot(3,1,1);myplot(B,A);%调用绘图函数myplot绘制损耗函数曲线yt='y_1(t)';subplot(3,1,2); tplot(y1t,T,yt);%调用绘图函数tplot绘制滤波器输出波形%===========================================%带通滤波器设计与实现=========================fpl=440; fpu=560; fsl=275; fsu=900;wp=[2*fpl/Fs,2*fpu/Fs];ws=[2*fsl/Fs,2*fsu/Fs]; rp=0.1; rs=60;[N,wp]=ellipord(wp,ws,rp,rs);%调用ellipord计算椭圆DF阶数N和通带截止频率wp[B,A]=ellip(N,rp,rs,wp);%调用ellip计算椭圆带通DF系统函数系数向量B和Ay2t=filter(B,A,st); %滤波器软件实现% 带通滤波器设计与实现绘图部分（省略）%补充1figure(2); subplot(3,1,3);myplot(B,A);%--figure(3);yt='y_2(t)';subplot(3,1,1); tplot(y2t,T,yt);%===========================================%高通滤波器设计与实现===========================fp=890; fs=600;wp=2*fp/Fs; ws=2*fs/Fs; rp=0.1; rs=60;%DF指标（低通滤波器的通、阻带边界频率）[N,wp]=ellipord(wp,ws,rp,rs);%调用ellipord计算椭圆DF阶数N和通带%截止频率wp [B,A]=ellip(N,rp,rs,wp,'high');%调用ellip计算椭圆带通DF系统函数系数向量B和Ay3t=filter(B,A,st); %滤波器软件实现% 高低通滤波器设计与实现绘图部分（省略）%补充2figure(3); subplot(3,1,2);myplot(B,A);%--figure(3);yt='y_3(t)';subplot(3,1,3); tplot(y3t,T,yt);%===========================================function st=mstg%产生信号序列向量st,并显示st的时域波形和频谱%st=mstg 返回三路调幅信号相加形成的混合信号,长度N=1600 N=1600 %N为信号st的长度Fs=10000;T=1/Fs;Tp=N*T;%采样频率Fs=10 kHz, Tp为采样时间t=0:T:(N-1)*T;k=0:N-1;f=k/Tp;fc1=Fs/10;%第1路调幅信号的载波频率fc1=1000 Hzfm1=fc1/10;%第1路调幅信号的调制信号频率fm1=100 Hzfc2=Fs/20;%第2路调幅信号的载波频率fc2=500 Hzfm2=fc2/10;%第2路调幅信号的调制信号频率fm2=50 Hzfc3=Fs/40;%第3路调幅信号的载波频率fc3=250 Hzfm3=fc3/10;%第3路调幅信号的调制信号频率fm3=25 Hzxt1=cos(2*pi*fm1*t).*cos(2*pi*fc1*t);%产生第1路调幅信号xt2=cos(2*pi*fm2*t).*cos(2*pi*fc2*t);%产生第2路调幅信号xt3=cos(2*pi*fm3*t).*cos(2*pi*fc3*t);%产生第3路调幅信号st=xt1+xt2+xt3; %三路调幅信号相加fxt=fft(st,N); %计算信号st的频谱%以下为绘图部分，绘制st的时域波形和幅频特性曲线subplot(3,1,1)plot(t,st);grid;xlabel('t/s');ylabel('s(t)');axis([0,Tp/8,min(st),max(st)]);title('(a) s(t)的波形')subplot(3,1,2)stem(f,abs(fxt)/max(abs(fxt)),'.');grid;title('(b) s(t)的频谱')axis([0,Fs/5,0,1.2]);xlabel('f/Hz');ylabel('幅度')function myplot(B,A)%时域离散系统损耗函数绘图%B为系统函数分子多项式系数向量%A为系统函数分母多项式系数向量[H,W]=freqz(B,A,1000);m = abs(H);plot(W/pi,20*log10(m/max(m))); grid on;xlabel('\omega/\pi'); ylabel('幅度(dB)');axis([0,1,-80,5]); title('损耗函数曲线');function tplot(xn,T,yn)%时域序列连续曲线绘图函数%xn：信号数据序列，yn：绘图信号的纵坐标名字(字符串) %T为采样间隔n = 0:length(xn)-1; t=n*T;plot(t,xn);xlabel('t/s'); ylabel(yn);axis([0,t(end),min(xn),1.2*max(xn)])实验五%实验五程序exp5.m%FIR数字滤波器设计及软件实现clear all; close all;%==调用xtg产生信号xt, xt长度N=1000,并显示xt及其频谱,=======N=1000;xt=xtg(N);fp=120; fs=150; Rp=0.2; As=60; Fs=1000;%输入给定指标%(1) 用窗函数法设计滤波器wc=(fp+fs)/Fs; %理想低通滤波器截止频率(关于pi归一化）B=2*pi*(fs-fp)/Fs; %过渡带宽度指标Nb=ceil(11*pi/B); %blackman窗的长度Nhn=fir1(Nb-1,wc,blackman(Nb));Hw=abs(fft(hn,1024));%求设计的滤波器频率特性ywt=fftfilt(hn,xt,N); %调用函数fftfilt对xt滤波%以下为用窗函数法设计法的绘图部分（滤波器损耗函数, 滤波器输出信号波形）（省略）f=[0:1023]*Fs/1024;figure(2)subplot(2,1,1)plot(f,20*log10(Hw/max(Hw)));grid;title('(3)低通滤波器幅频特性')axis([0,Fs/2,-120,20]);xlabel('f/Hz');ylabel('幅度')t=[0:2*N-1]/Fs;Tp=N/Fs; %%%%%%%%%%subplot(2,1,2);plot(t,ywt);grid;axis([0,Tp/2,-1,1]);xlabel('t/s');ylabel('y_w(t)');title('(4)滤波噪声后的信号波形');%(2) 用等波纹最佳逼近法设计滤波器fb=[fp,fs]; m=[1,0];%确定remezord函数所需参数f,m,devdev=[(10^(Rp/20)-1)/(10^(Rp/20)+1),10^(-As/20)];[Ne,fo,mo,W]=remezord(fb,m,dev,Fs);%确定remez函数所需参数hn=remez(Ne,fo,mo,W);%调用remez函数进行设计Hw=abs(fft(hn,1024));%求设计的滤波器频率特性yet=fftfilt(hn,xt,N); %调用函数fftfilt对xt滤波%以下为用等波纹设计法的绘图部分（滤波器损耗函数, 滤波器输出信号yw(t)波形）%(省略)f=[0:1023]*Fs/1024;figure(3)subplot(2,1,1)plot(f,20*log10(Hw/max(Hw)));grid;title('(5)低通滤波器幅频特性')axis([0,Fs/2,-80,10]);xlabel('f/Hz');ylabel('幅度')%t=[0:N-1]/Fs;Tp=N/Fs;subplot(2,1,2)plot(t,yet);grid;axis([0,Tp/2,-1,1]);xlabel('t/s');ylabel('y_e(t)');title('(6)滤波噪声后的信号波形');function xt=xtg(N)%实验五信号x(t)产生函数,并显示信号的幅频特性曲线%xt=xtg 产生一个长度为N,有加性高频噪声的单频调幅信号xt,采样频率Fs=1000 Hz %载波频率fc=Fs/10=100Hz，调制正弦波频率f0=fc/10=10 HzN=2000; Fs=1000; T=1/Fs; Tp=N*T;t=0:T:(N-1)*T;fc=Fs/10; f0=fc/10; %载波频率fc=Fs/10, 单频调制信号频率为f0=Fc/10mt=cos(2*pi*f0*t);%产生单频正弦波调制信号mt, 频率为f0ct=cos(2*pi*fc*t); %产生载波正弦波信号ct, 频率为fcxt=mt.*ct; %相乘产生单频调制信号xtnt=2*rand(1,N)-1; %产生随机噪声nt%====设计高通滤波器hn,用于滤除噪声nt中的低频成分,生成高通噪声====fp=150; fs=200; Rp=0.1; As=70;% 滤波器指标fb=[fp,fs]; m=[0,1];% 计算remezord函数所需参数f,m,devdev=[10^(-As/20),(10^(Rp/20)-1)/(10^(Rp/20)+1)];[n,fo,mo,W]=remezord(fb,m,dev,Fs);% 确定remez函数所需参数hn=remez(n,fo,mo,W);% 调用remez函数进行设计，用于滤除噪声nt中的低频成yt=filter(hn,1,10*nt);%滤除随机噪声中低频成分, 生成高通噪声yt%=======================================xt=xt+yt; %噪声加信号fst=fft(xt,N); k=0:N-1; f=k/Tp;subplot(3,1,1); plot(t,xt); grid;xlabel('t/s'); ylabel('x(t)');axis([0,Tp/5,min(xt),max(xt)]);title('(a) 信号加噪声波形')subplot(3,1,2); plot(f,abs(fst)/max(abs(fst)));grid; title('(b) 信号加噪声的频谱')axis([0,Fs/2,0,1.2]);xlabel('f/Hz');ylabel('幅度');实验6%实验六程序exp6.m% DTMF双频拨号信号的生成和检测程序%clear all;clc;tm=[1,2,3,65;4,5,6,66;7,8,9,67;42,0,35,68]; % DTMF信号代表的16个数N=205;K=[18, 20, 22, 24, 31, 34, 38, 42];f1=[697,770,852,941]; % 行频率向量f2=[1209,1336,1477,1633]; % 列频率向量TN=input('键入6位电话号码= '); % 输入6位数字TNr=0; % 接收端电话号码初值为零for l=1:6;d=fix(TN/10^(6-l));TN=TN-d*10^(6-l);for p=1:4;for q=1:4;if tm(p,q)==abs(d); break,end% 检测码相符的列号qendif tm(p,q)==abs(d); break,end% 检测码相符的行号pendn=0:1023; % 为了发声, 加长序列x = sin(2*pi*n*f1(p)/8000) + sin(2*pi*n*f2(q)/8000);% 构成双频信号sound(x,8000); % 发出声音pause(0.1)% 接收检测端的程序X=goertzel(x(1:205),K+1);% 用Goertzel算法计算八点DFT样本val = abs(X); % 列出八点DFT向量subplot(3,2,l);stem(K,val,'.');grid;xlabel('k');ylabel('|X(k)|') % 画出DFT(k)幅度axis([10 50 0 120])limit = 80;for s=5:8;if val(s) > limit, break, end % 查找列号endfor r=1:4;if val(r) > limit, break, end % 查找行号endTNr=TNr+tm(r,s-4)*10^(6-l);enddisp('接收端检测到的号码为:')% 显示接收到的字符disp(TNr)%实验六程序exp6.m%改成识别8位电话号码% DTMF双频拨号信号的生成和检测程序%clear all;clc;tm=[1,2,3,65;4,5,6,66;7,8,9,67;42,0,35,68]; % DTMF信号代表的16个数N=205;K=[18, 20, 22, 24, 31, 34, 38, 42];f1=[697,770,852,941]; % 行频率向量f2=[1209,1336,1477,1633]; % 列频率向量TN=input('键入8位电话号码= '); % 输入8位数字TNr=0; % 接收端电话号码初值为零for l=1:8;d=fix(TN/10^(8-l));TN=TN-d*10^(8-l);for p=1:4;for q=1:4;if tm(p,q)==abs(d); break,end% 检测码相符的列号qendif tm(p,q)==abs(d); break,end% 检测码相符的行号pendn=0:1023; % 为了发声, 加长序列x = sin(2*pi*n*f1(p)/8000) + sin(2*pi*n*f2(q)/8000);% 构成双频信号sound(x,8000); % 发出声音pause(0.1)% 接收检测端的程序X=goertzel(x(1:205),K+1);% 用Goertzel算法计算八点DFT样本val = abs(X); % 列出八点DFT向量subplot(4,2,l);stem(K,val,'.');grid;xlabel('k');ylabel('|X(k)|') % 画出DFT(k)幅度axis([10 50 0 120])limit = 80;for s=5:8;if val(s) > limit, break, end % 查找列号endfor r=1:4;if val(r) > limit, break, end % 查找行号endTNr=TNr+tm(r,s-4)*10^(8-l);enddisp('接收端检测到的号码为:')% 显示接收到的字符disp(TNr)。