《数字信号处理上机实验代码..pdf》由会员分享,可在线阅读,更多相关《数字信号处理上机实验代码..pdf(19页珍藏版)》请在三一文库上搜索。
1、文件名 :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
2、; plot(t,xn; xlabel(t/s;ylabel(yn; axis(0,t(end,min(xn,1.2*max(xn; 文件名 :myplot.m (实验一、四需要 %(1myplot;计算时域离散系统损耗函数并绘制曲线图。 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(损耗函数曲线 ;
3、文件名 :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 (实验一需要 %(2mpplot;计算时域离散系统损耗函数和相频特性函数, 并绘制曲 线图。 funct
4、ion mpplot(B,A,Rs %mpplot(B,A,Rs %时域离散系统损耗函数和相频特性绘图 %B 为系统函数分子多项式系数向量 %A 为系统函数分母多项式系数向量 %Rs为滤波器阻带最小衰减 ,省略则幅频曲线最小值取-80dB if nargin3 ymin=- 80;else ymin=-Rs-20;end;% 确定幅频曲线 纵坐标最小值 H,W=H,W=freqz(B,A,1000 m=abs(H; subplot(2,2,1; plot(W/pi,20*log10(m/max(m;grid on; xlabel(omega/pi;ylabel(幅度 (dB axis(0,1,
5、ymin,5;title( 损耗函数曲线 ; subplot(2,2,3; plot(W/pi,p/pi; xlabel(omega/pi;ylabel(相位 /pi;grid on; title(b 相频特性曲线 ; 文件名 :mfftplot.m ( 实验一需要 function mfftplot(xn,N %mfftplot(xn,N 计算序列向量 xn 的 N 点 fft 并绘制其幅频特性曲线 Xk=fft(xn,N;% 计算信号 xn 的频谱的 N 点采样 %=以下为绘图部分 = k=0:N-1;wk=2*k/N; m=abs(Xk;mm=max(m; plot(wk,m/mm;gr
6、id on; xlabel(omega/pi;ylabel(幅度 (dB; axis(0,2,0,1.2; title(幅度特性曲线 ; 文件名 :mstg.m (实验四需要 程序: function st=mstg %产生信号序列向量 st, 并显示 st 的时域波形和频谱 %st=mstg 返回三路调幅信号相加形成的混合信号,长度 N=1600 N=1600; %N 为 信号 st 的长度 Fs=10000;T=1/Fs;Tp=N*T; %采样频率 Fs=10kHz, Tp 为采样 时间 t=0:T:(N-1*T;k=0:N-1;f=k/Tp; fc1=Fs/10; %第 1 路调幅信号的
7、载波频率 fc1=1000Hz fm1=fc1/10; %第 1路调幅信号的调制信号频率 fm1=100Hz fc2=Fs/20; %第 2 路调幅信号的载波频率 fc2=500Hz fm2=fc2/10; %第 2路调幅信号的调制信号频率 fm2=50Hz fc3=Fs/40; %第 3路 调幅信号的载波频率 fc3=250Hz fm3=fc3/10; %第 3 路调幅信号的调制信号频率 fm3=25Hz xt1=cos(2*pi*fm1*t.*cos(2*pi*fc1*t; %产生第 1 路调幅信号 xt2=cos(2*pi*fm2*t.*cos(2*pi*fc2*t; %产生第 2路调幅
8、信号 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 的频谱
9、 ; axis(0,Fs/5,0,1.2; xlabel(f/Hz;ylabel( 幅度 ; 文件名 :xtg.m (实验五需要 程序: function xt=xtg(N %实验五信号 x(t 产生 , 并显示信号的幅频特性曲线 %xt=xtg(N 产生一个长度为 N, 有加性高频噪声的单频调幅信号 xt, 采样频率 Fs=1000Hz %载波频率 fc=Fs/10=100Hz,调制正弦波频率 f0=fc/10=10Hz. Fs=1000;T=1/Fs;Tp=N*T; t=0:T:(N-1*T; fc=Fs/10;f0=fc/10; %载波频率 fc=Fs/10,单频调制信号频率为 f0=F
10、c/10; mt=cos(2*pi*f0*t; % 产生单频正弦波调制信号 mt ,频率为 f0 ct=cos(2*pi*fc*t; % 产生载波正弦波信号 ct ,频率为 fc xt=mt.*ct; % 相乘产生单频调制信号 xt nt=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,dev dev=10(-As/20,(10(Rp/20-1/(10(
11、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 信
12、号加噪声波形 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(幅度 10.1 系统响应及系统稳定性 close all;clear all;clc; %内容 1:调用 filter 解差分方程 ,由系统对 u(n的响应判断稳定 性 %= A=1,-0.9;B=0.05,0.05; % 系统差分方程系数向量 B 和 A x1n=1 1 1 1 1 1 1 1 zeros(1,50; %产生信号 x1(n=R8(n x2n=ones(1,128;
13、 % 产生信号 x2(n=u(n hn=impz(B,A,58; %求系统单位脉冲响应 h(n subplot(2,2,1;y=h(n;tstem(hn,y; % 调用函数 tstem 绘图 title(a 系统单位脉冲响应 h(n;box on y1n=filter(B,A,x1n; % 求系统对 x1(n 的响应 y1(n subplot(2,2,2;y=y1(n;tstem(y1n,y; title(b 系统对 R8(n的响应 y1(n;box on y2n=filter(B,A,x2n; % 求系统对 x2(n 的响应 y2(n subplot(2,2,4;y=y2(n;tstem(y
14、2n,y; title(c 系统对 u(n的响应 y2(n;box on %内容 2:调用 conv 函数计算卷积 %= x1n=1 1 1 1 1 1 1 1 ; %产生信号 x1(n=R8(n h1n=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;box on subplot(2,2,2;y=y21(n
15、;tstem(y21n,y; title(e h1(n 与 R8(n的卷积 y21(n;box on subplot(2,2,3;y=h2(n;tstem(h2n,y; % 调 用函数 tstem 绘图 title(f 系统单位脉冲响应 h2(n;box on subplot(2,2,4;y=y22(n;tstem(y22n,y; title(g h2(n 与 R8(n的卷积 y22(n;box on %内容 3:谐振器分析 %= un=ones(1,256; % 产生信号 u(n n=0:255; xsin=sin(0.014*n+sin(0.4*n; %产生正弦信号 A=1,-1.8237
16、,0.9801;B=1/100.49,0,-1/100.49; %系统差分方 程系数向量 B 和 A y31n=filter(B,A,un; % 谐振器对 u(n的响应 y31(n y32n=filter(B,A,xsin; % 谐振器对 u(n的响应 y31(n figure(3 subplot(2,1,1;y=y31(n;tstem(y31n,y; title(h 谐振器对 u(n 的响应 y31(n;box on subplot(2,1,2;y=y32(n;tstem(y32n,y; title(i 谐振器对正弦信号的响应 y32(n;box on 10.2 时域采样与频域采样 时域采样
17、理论验证程序 exp2a.m : %= close all;clear all;clc; Tp=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*20.5;omega=pi*50*20.5; xnt=A*exp(-alph*n*T.*sin(omega*n*T; Xk=T*fft(xnt,M; %M点 FFTxnt yn=xa(nT;subplot(3,2,1; tstem(xnt,yn; %调用自编绘图函数 tstem 绘
18、制序列图 box on;title(a Fs=1000Hz; k=0:M-1;fk=k/Tp; subplot(3,2,2;plot(fk,abs(Xk;title(a T*FTxa(nT,Fs=1000Hz; xlabel(f(Hz;ylabel( 幅度 ;axis(0,Fs,0,1.2*max(abs(Xk %= % Fs=300Hz和 Fs=200Hz的程序与上面 Fs=1000Hz完全相同。 频域采样理论 验证程序 exp2b.m : %= close all;clear all;clc; M=27;N=32;n=0:M; %产生 M 长三角波序列 x(n xa=0:floor(M/2
19、; xb= ceil(M/2-1:-1:0; xn=xa,xb; Xk=fft(xn,1024; %1024 点 FFTx(n, 用于近似序列 x(n 的 TF X32k=fft(xn,32 ; %32 点 FFTx(n x32n=ifft(X32k; %32 点 IFFTX32(k 得到 x32(n X16k=X32k(1:2:N; % 隔点抽取 X32k 得到 X16(K x16n=ifft(X16k,N/2; %16 点 IFFTX16(k 得到 x16(n subplot(3,2,2;stem(n,xn,.;box on title(b 三角波序列 x(n;xlabel(n;ylabe
20、l(x(n;axis(0,32,0,20 k=0:1023;wk=2*k/1024; subplot(3,2,1;plot(wk,abs(Xk;title(aFTx(n; xlabel(omega/pi;ylabel(|X(ejomega|;axis(0,1,0,200 k=0:N/2-1; subplot(3,2,3;stem(k,abs(X16k,.;box on title(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 on title(d
21、 16 点 IDFTX_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 on title(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 on title(f 32 点 IDFTX_3_2(k;xlabel(n;ylabel(x_3_2(n;axis(0,32,0,20 10.3 用 FFT 对信号做 频谱分析 exp3
22、01.m : %= clear all;close all;clc; %实验内容 (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点 DFT X1k16=fft(x1n,16; % 计算 x1n 的 16点 DFT X2k8=fft(x2n,8; % 计算 x1n 的 8点 DFT X2k16=fft(x2n,16; % 计算 x1n 的 16点 DFT X3k8=ff
23、t(x3n,8; % 计算 x1n 的 8点 DFT X3k16=fft(x3n,16; % 计算 x1n 的 16点 DFT %以下绘制幅频特性曲线 subplot(2,2,1;mstem(X1k8; %绘制 8 点 DFT 的幅频特性图 title(1a 8 点 DFTx_1(n;x label(/ ;ylabel(幅度 ; axis(0,2,0,1.2*max(abs(X1k8 subplot(2,2,3;mstem(X1k16; %绘制 16点 DFT 的幅频特性图 title(1b16 点 DFTx_1(n;xlabel(/ ;ylabel(幅度 ; axis(0,2,0,1.2*m
24、ax(abs(X1k16 figure(2 subplot(2,2,1;mstem(X2k8; %绘制 8 点 DFT 的幅频特性图 title(2a 8 点 DFTx_2(n;xlabel(/ ;ylabel(幅度 ; axis(0,2,0,1.2*max(abs(X2k8 subplot(2,2,2;mstem(X2k16; %绘制 16点 DFT 的幅频特性图 title(2b16 点 DFTx_2(n;xlabel(/ ;ylabel(幅度 ; axis(0,2,0,1.2*max(abs(X2k16 subplot(2,2,3;mstem(X3k8; %绘制 8点 DFT 的幅频特性
25、图 title(3a 8 点 DFTx_3(n;xlabel(/ ;ylabel(幅度 ; axis(0,2,0,1.2*max(abs(X3k8 subplot(2,2,4;mstem(X3k16;%绘制 16点 DFT 的幅频特性图 title(3b16 点 DFTx_3(n;xlabel(/ ;ylabel(幅度 ; exp302.m : %实验内容 (2 周期序列谱分析 %= close all;clear all;clc; N=8;n=0:N-1; %FFT 的变换区间 N=8 x4n=cos(pi*n/4; x5n=cos(pi*n/4+cos(pi*n/8; X4k8=fft(x
26、4n; % 计算 x4n 的 8点 DFT X5k8=fft(x5n; % 计算 x5n 的 8点 DFT N=16;n=0:N-1; %FFT 的变换区间 N=16 x4n=cos(pi*n/4; x5n=cos(pi*n/4+cos(pi*n/8; X4k16=fft(x4n; % 计算 x4n 的 16点 DFT X5k16=fft(x5n; % 计算 x5n 的 16点 DFT figure(3 subplot(2,2,1;mstem(X4k8; %绘制 8 点 DFT 的幅频特性图 title(4a 8 点 DFTx_4(n;xlabel(/ ;ylabel(幅度 ; axis(0,
27、2,0,1.2*max(abs(X4k8 subplot(2,2,3;mstem(X4k16; %绘制 16点 DFT 的幅频特性图 title(4b16 点 DFTx_4(n;xlabel(/ ;ylabel(幅度 ; subplot(2,2,2;mstem(X5k8; %绘制 8 点 DFT 的幅频特性图 title(5a 8 点 DFTx_5(n;xlabel(/ ;ylabel(幅度 ; axis(0,2,0,1.2*max(abs(X5k8 subplot(2,2,4;mstem(X5k16; %绘制 16点 DFT 的 幅频特性图 title(5b16 点 DFTx_5(n;xla
28、bel(/ ;ylabel(幅度 ; axis(0,2,0,1.2*max(abs(X5k16 exp303.m : %实验内容 (3 模拟周期信号谱分析 %= close all;clear all;clc; figure(4 Fs=64;T=1/Fs; N=16;n=0:N-1; %FFT 的变换区间 N=16 x6nT=cos(8*pi*n*T+cos(16*pi*n*T+cos(20*pi*n*T; %对 x6(t16 点采样 X6k16=fft(x6nT; % 计算 x6nT 的 16 点 DFT X6k16=fftshift(X6k16; % 将零频率移到频谱中心 Tp=N*T;F
29、=1/Tp; % 频率分辨率 F k=-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点 |DFTx_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=16 x6nT=cos(8*pi*n*T+cos(16*pi*n*T+cos(20*pi*n*T
30、; %对 x6(t32 点采样 X6k32=fft(x6nT; % 计算 x6nT 的 32 点 DFT X6k32=fftshift(X6k32; % 将零频率移到频谱中心 Tp=N*T;F=1/Tp; % 频率分辨率 F k=-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 点 |DFTx_6(nT|;xlabel(f(Hz;ylabel( 幅度 ; axis(-N*F/2-1,N*F/2- 1
31、,0,1.2*max(abs(X6k32 N=64;n=0:N-1; %FFT 的变换区间 N=16 x6nT=cos(8*pi*n*T+cos(16*pi*n*T+cos(20*pi*n*T; %对 x6(t64 点采样 X6k64=fft(x6nT; % 计算 x6nT 的 64 点 DFT X6k64=fftshift(X6k64; % 将零频率移到频谱中心 Tp=N*T;F=1/Tp; % 频率分辨率 F k=-N/2:N/2-1;fk=k*F; % 产生 16点 DFT 对应的采样点频率 (以 零频率为中心 subplot(3,1,3;stem(fk,abs(X6k64,.; box
32、 on% 绘制 8 点 DFT 的 幅频特性图 title(6a 64点 |DFTx_6(nT|;xlabel(f(Hz;ylabel( 幅度 ; axis(-N*F/2-1,N*F/2- 1,0,1.2*max(abs(X6k64 10.4 IIR 数字滤波器设计及软件实现 exp4.m : %= clear all;close all;clc; Fs=10000;T=1/Fs; % 采样频率 %调用信号产生函数 mstg 产生由三路抑制载波调幅信号相加构成的复合信号 st st=mstg; %低通滤波器设计与实现 %= fp=280;fs=450; wp=2*fp/Fs;ws=2*fs/F
33、s;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 和 A y1t=filter(B,A,st; % 滤波器软件实现 % 低通滤波器设计与实现绘图部分 figure(2;subplot(3,1,1; myplot(B,A; % 调用绘图函数 myplot 绘制损耗函数曲线 yt=y_1(t; subplot(3,1,2;tplot(
34、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 和 A y2t=filter(B,A,st; % 滤波器软件实现 % 带通滤波器设
35、计与实现绘图部分(省略 %高通滤波器设计与实现 %= 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 和 A y3t=filter(B,A,st; % 滤波器软件实现 % 高低通滤波器设计与实现绘图部分(省略 10.5 FIR数字滤波器设计及软件实现
36、 exp5.m : %= clear all;close all;clc; %调用 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 窗的长度 N hn=fir1(Nb-1,wc,blackman(Nb; Hw=abs(ff
37、t(hn,1024; % 求设计的滤波器频率特性 ywt=fftfilt(hn,xt,N; % 调用函数 fftfilt 对 xt 滤波 %以下为用窗函数法设计法的绘图部分(滤波器损耗函数 ,滤波器输 出信号波形 %省略 % (2 用等波纹最佳逼近法设计滤波器 %= fb=fp,fs;m=1,0; % 确定 remezord 函数所需参数 f,m,dev dev=(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(nT 波形) %省略
链接地址:https://www.31doc.com/p-4954667.html