数字信号实验

实验程序：
1 产生10点的单位抽样序列δ(n)；
%function unit_pulse(N)
% unit_pulse.m
N=10；
x=zeros(1,N);
x(1)=1;
n=0:N-1;
figure(1);
stem(n,x);xlabel('单位抽样序列')
axis([-1 20 0 1.1])
2产生10点同时移位3位的单位抽样序列δ(n-3)；
%function shift_unit_pulse(N,k)
% shift_unit_pulse.m
N=10;
k=3;
x=zeros(1,N);
x(k+1)=1;
n=0:N-1;
figure(2);
stem(n,x);xlabel('移位3位的单位抽样序列')
axis([-1 20 0 1.1])
或function [x, n]=i shift_unit_pulse (n0,ns,nf)
n=[0:9]；
x=[(n-3)==0]
3产生任意序列
f(n)=8δ(n)+7δ(n-1)+6δ(n-2) +5δ(n-3)+ 4δ(n-4)+7δ(n-5)；%function arbitrary_pulse(N)
% arbitrary_pulse.m
N=10
x=zeros(1,N);
x(1)=8;x(2)=7;x(3)=6;x(4)=5;x(5)=4;x(6)=7;
n=0:N-1;
figure(3);
stem(n,x);xlabel('任意序列f(n)')
axis([-1 20 0 9])
4产生N＝32点的单位阶跃序列
%function unit_step(N)
% unit_step.m
N=32;
x=ones(1,N);
n=0:N-1;
figure(4);
stem(n,x);xlabel('单位阶跃序列') axis([-1 32 0 1.1])
5产生斜率为3，n
0=4,点数为20点的斜坡序列g(n)=B(n-n
)
%function slope(N,k,B)
% slope.m
N=20;k=4;B=3;
x=[zeros(1,k) ones(1,N-k)];
for i=1:N
x(i)=B*x(i)*(i-k);
end
n=0:N-1;
figure(5);
stem(n,x);xlabel('斜坡序列')
axis([-1 10 0 90])
6产生幅度A=3，频率f＝100，初始相位 ＝1.2，点数为32点的正弦序列。

%function sine(N,A,f,fai)
% sine.m
A=3;f=100;fai=1.2;N=32;
n=0:N-1;
x=A*sin(2*pi*f*(n/N)+fai);
figure(6);
stem(n,x);xlabel('正弦序列')
axis([-1 32 -3.2 3.2])
7产生幅度A=3，角频率ω＝314，点数为32点的复正弦序列。

%function complex_sine(N,A,w)
% complex_sine.m
A=3; w=314;N=32
n=0:N-1;
x=A*exp(j*w*n);
figure(7);
stem(n,x);xlabel('复正弦序列')
axis([-1 32 -3.2 3.2])
8产生幅度A=3，a＝0.7，点数为32点的实指数序列。

%function real_exponent(N,A,a)
% real_exponent.m
a=0.7;A=3;
n=0:N-1;
x=A*a.^n;
figure(8);
stem(n,x);xlabel('实指数序列')
axis([-1 32 0 3.2])
四：实验结果和分析
单位抽样序列
移位3位的单位抽样序列
任意序列f(n)
单位阶跃序列
斜坡序列
正弦序列
复正弦序列
实指数序列
实验程序： 1
求序列()cos(0.225)x n n π=的谱分析。

160长度更准确，整周期采样估计更准确。

clear all close all
clc
%序列周期为r(2*pi/w)=r(2/0.225)=80/9 %选定序列的长度为80点可准确估计 N=160;
n=0:N-1;
x=cos(0.225*pi*n);%生成信号 y=fft(x); %求信号的频谱 subplot(2,1,1) %画图 plot(n,x) xlabel('n') ylabel('x(n)')
subplot(2,1,2) plot(n/N*2,abs(y))
xlabel('nomorlized frequency') ylabel('Y(exp(j\omega))') xlim([0,1]) figure;
subplot(2,1,1) stem(n,x) xlabel('n') ylabel('x(n)') subplot(2,1,2) stem(n,abs(y)) xlabel('k') ylabel('Y(k))')
N=256; %点数不同，估计准确度不同 n=0:N-1;
x=cos(0.225*pi*n); y=fft(x); figure; subplot(2,1,1) plot(n,x) xlabel('n') ylabel('x(n)') subplot(2,1,2)
plot(n/N*2,abs(y))
xlabel('nomorlized frequency') ylabel('Y(exp(j\omega))')
xlim([0,1]) figure;
subplot(2,1,1) stem(n,x) xlabel('n') ylabel('x(n)') subplot(2,1,2) stem(n,abs(y)) xlabel('k') ylabel('Y(k))')
2求信号()cos()cos(2.5)x t t t ππ=+的谱分析，绘出()X j Ω的幅相特性。

给出采样频率s f 、采样点数N ，谱分辨率F 等参数的选择原因及结果。

clear all close all clc
f1=0.5; f2=1.25;
fs=12.5;ts=1/fs;F=0.25;Tp=1/F; t=0:ts:Tp;
xt=cos(2*pi*f1*t)+cos(2*pi*f2*t); y=fft(xt); %求信号的频谱 plot(t,xt)
xlabel('time/s') ylabel('x(t)') title('signal') figure
subplot(2,1,1) %画图 plot((0:length(y)-1)*F,abs(y)) xlabel('frequency/Hz') ylabel('Y(j\Omega)')
subplot(2,1,2)
plot((0:length(y)-1)*F,angle(y)) xlabel('frequency/Hz') ylabel('\Phi(j\Omega)')
尝试变换参数，分析信号频谱。

3 观察并分析采用不同频率时，对函数50()218.2sin(50)()t a x t e t u t ππ-=的频谱影响。

（a ）：以1000s f Hz =，对其进行采样得到x 1(n)。

（b ）：以200s f Hz =，对其进行采样得到x 3(n) a=218.2;
b=50*pi;
fs=1000; %采样频率1000hz T=1/fs;
Tp=50*0.001; %观察时间50微秒
M=Tp*fs; %采样点数 n=0:M-1
y=a*exp(-b*n*T).*sin(b*n*T ) %函数表达式 subplot(3,2,1) stem(n,y)
xlabel('n');ylabel('xa(nT)');title('fs=1000HZ'); axis([0,M,-2,1.2*max(abs(y))])
yk=T*fft(y,M) %M 点FFT K=0:M-1; fk=K/Tp;
subplot(3,2,2) plot(fk,abs(yk))
xlabel('f(Hz)');ylabel('幅度');title('T*FFT,fs=1000HZ'); axis([0,fs,0,1.2*max(abs(yk))])
fs=200; %采样频率200hz T=1/fs;
Tp=50*0.001; M=Tp*fs; n=0:M-1
y=a*exp(-b*n*T).*sin(b*n*T )
subplot(3,2,5) stem(n,y)
xlabel('n');ylabel('xa(nT)');title('fs=200HZ'); axis([0,M,-2,1.2*max(abs(y))]);
yk=T*fft(y,M) %M 点FFT K=0:M-1;
fk=K/Tp; subplot(3,2,6) plot(fk,abs(yk))
xlabel('f(Hz)');ylabel('幅度');title('T*FFT,fs=200HZ'); axis([0,fs,0,1.2*max(abs(yk))])
四：实验结果和分析
n
x (n
)
nomorlized frequency
Y (e x p (j ))
n
x (n )
k
Y (k ))
n
x (n
)
nomorlized frequency
Y (e x p (j ω))
n
x (n
)
k
Y (k ))
time/s
x (t )
frequency/Hz
Y (j Ω)
frequency/Hz
Φ(j Ω)
n x a (n T )
f(Hz)
幅度
n
x a (n T )
f(Hz)
幅度
三、实验所采用的功能函数
1.巴特沃斯滤波器阶数选择函数
（1）[N,wc]=buttord(wp，ws，αp，αs)
作用：计算巴特沃斯数字滤波器的阶数N和3dB截止频率wc, wc为数字频率，单位rad。

说明：调用参数wp，ws分别为数字滤波器的通带、阻带截止频率的归一化值，要求：0≤wp≤1，0≤ws≤1。

αp，αs分别为通带最大衰减和组带最小衰减（dB）。

当ws≤wp时，为高通滤波器；当wp和ws为二元矢量时，为带通或带阻滤波器，这时wc也是二元向量。

（2）[N,Ωc]=buttord(Ωp，Ωs，αp，αs，‘s’)
作用：计算巴特沃斯模拟滤波器的阶数N和3dB截止频率Ωc。

说明：Ωp，Ωs，Ωc均为实际模拟角频率。

模拟频率f：每秒经历多少个周期，单位Hz，即1/s，信号的真实频率，可用于模拟信号和数字信号；
模拟角频率Ω：每秒经历多少弧度，单位rad/s，通常只于模拟信号；
数字频率w：每个采样点间隔之间的弧度，单位rad，通常只用于数字信号。

关系：Ω=2pi*f；w = Ω*T=2pi*f/F。

(F=1/Ts为采样频率，Ts为采样间隔)
2.完整巴特沃斯滤波器设计函数
（1）格式： [b，a]=butter（N，wc，‘ftype’）
作用：计算N阶巴特沃斯数字滤波器系统函数分子、分母多项式的系数向量b、a。

说明：调用参数N和wc分别为巴特沃斯数字滤波器的阶数和3dB截止频率的归一化值，一般是调用buttord格式（1）计算N和wc。

系数b、a是按照z-1的升幂排列。

（2）格式：[B，A]=butter（N，Ωc，‘ftype’，‘s’）
作用：计算巴特沃斯模拟滤波器系统函数的分子、分母多项式系数向量。

说明：调用参数N和Ωc分别为巴特沃斯模拟滤波器的阶数和3dB截止频率(实
际角频率)，可调用buttord （2）格式计算N 和Ωc 。

系数B 、A 按s 的正降幂排列。

tfype 为滤波器的类型：
◇ftype=high 时，高通；Ωc 只有1个值。

◇ftype=stop 时，带阻；Ωc=[Ωcl,Ωcu]，分别为带阻滤波器的通带3dB 下截止频率和上截止频率。

◇ ftype 缺省时：若Ωc 只有1个值，则默认为低通；若Ωc 有2个值，则默认为带通；其通带频率区间Ωcl<Ω<Ωcu 。

3. 求离散系统频响特性的函数freqz() 格式：[H,w]=freqz(b,a,N)
说明：b 和a 分别为离散系统的系统函数分子、分母多项式的系数向量，返回量H 则包含了离散系统频响在 0～pi 范围内N 个频率等分点的值（其中N 为正整数），w 则包含了范围内N 个频率等分点。

调用默认的N 时，其值是512。

可以先调用freqz()函数计算系统的频率响应，然后利用abs()和angle()函数及plot()函数，绘制出系统的频响曲线。

4. 滤波器离散化函数：
bilinear （使用双线性变换法把模拟滤波器转换为数字滤波器） impinvar （使用脉冲响应不变法把模拟滤波器转换为数字滤波器）
四、实验内容及步骤
1. 用直接设计法设计BW （巴特沃斯）低通数字滤波器。

采样频率为2000Hz ，通带中允许的最大衰减为0.5dB ，阻带内的最小衰
减为40dB ，通带上限临界频率为30Hz ，阻带下限临界频率为40Hz 。

设计步骤：
（1） 确定滤波器的设计指标：s p s ααωω、、、p ；
（2） 运用函数s p s buttord ,,,p ωωαα（）计算巴特沃斯低通滤波器的阶
数N 和归一化3db 截止频率c ω；
（3） 运用函数butter ,c N ω（）求得低通滤波器的系统函数的分子、分
母多项式形式；
（4） 作图显示滤波器的幅频特性和相位特性。

2. 脉冲响应不变法设计数字滤波器
使用脉冲响应不变法设计数字低通滤波器，其指标为：通带临界频率0.5π，通带内衰减小于1dB ；阻带临界频率0.8π，阻带内衰减大于15dB,采样频率为100Hz 。

设计步骤：
（1） 确定数字频率指标；
（2） 采用脉冲响应不变法求得模拟低通滤波器频率设计指标； （3） 用butterworth 设计方法求得模拟低通滤波器的截止频率和阶
数；
（4） 设计归一化模拟低通滤波器；
（5） 利用脉冲响应不变法把模拟滤波器转换为数字滤波器； （6） 画出幅度响应和相位响应图。

3. 应用双线性变换方法设计低通数字滤波器
数字低通滤波器的设计指标为：通带临界频率0.5π，通带内衰减小于1dB ；阻带临界频率0.8π，阻带内衰减大于15dB,采样频率为100Hz 。

设计步骤：
（1） 确定数字频率指标；
（2） 采用双线性变换法求得模拟低通滤波器频率设计指标； （3） 用butterworth 设计方法求得模拟低通滤波器的截止频率和阶
数；
（4） 设计归一化模拟低通滤波器；
（5） 利用双线性变换法把模拟滤波器转换为数字滤波器； （6） 画出幅度响应和相位响应图。

实验程序：
1直接设计法设计BW （巴特沃斯）低通数字滤波器
fp=40; %带通截止频率
fs=30; %阻通截止频率
ft=200; %采样频率
rp=0.5;
rs=40;
wp=fp/(ft/2); %利用Nyquist频率进行归一化
ws=fs/(ft/2);
[n,wc]=buttord(wp,ws,rp,rs); %求数字滤波器的最小阶数和截止频率
[b,a]=butter(n,wc); %设计低通数字滤波器系数b，a
[H,W]=freqz(b,a); %求系统频响特性，W为数字角频率，单位rad figure;
plot(W*ft/(2*pi),abs(H));grid; %绘出频率响应曲线
xlabel('频率/Hz');ylabel('幅值');
figure
plot(W*ft/(2*pi),angle(H));grid; %绘出频率响应曲线
xlabel('频率/Hz');ylabel('相位');
2脉冲响应不变法设计数字滤波器
clear all;
fs=100
wp=0.5*pi;
ws=0.8*pi;
rp=1;
rs=15;
Wp=wp*fs; %由数字角频率转换为模拟角频率（脉冲响应不变法）
Ws=ws*fs;
[n,wc]=buttord(Wp,Ws,rp,rs,'s'); %选择滤波器的最小阶数
[b,a]=butter(n,wc,'s');
[bz,az]=impinvar(b,a,fs); %脉冲相应不变法变换为数字滤波器
[H,W]= freqz(bz,az); %求解数字滤波器的频率响应
plot(W*fs/(2*pi),abs(H));
grid;
xlabel('频率/hz');ylabel('幅值/dB');
title('脉冲响应不变变换法')
figure;
plot(W/pi,20*log10(abs(H)));
grid;
xlabel('归一化频率');ylabel('幅值/dB');
title('脉冲响应不变变换法')
3双线性变换法设计数字滤波器
clear all;closeall;clc
fs=100;
wp=0.5*pi;
ws=0.8*pi;
rp=1; rs=15;
Wp=2*fs*tan(wp/2); %由数字角频率转换为模拟角频率（脉冲响应不变法） Ws=2*fs*tan(ws/2);
[n,wc]=buttord(Wp,Ws,rp,rs,'s'); %选择滤波器的最小阶数
[b,a]=butter(n,wc,'s');
[bz,az]=bilinear(b,a,fs); %脉冲相应不变法变换为数字滤波器 [H,W]= freqz(bz,az); %求解数字滤波器的频率响应 plot(W*fs/(2*pi),20*log10(abs(H))); grid;
xlabel('频率/hz');ylabel('幅值/dB'); title('双线性变换法') figure; plot(W/pi,20*log10(abs(H))); grid;
xlabel('归一化频率');ylabel('幅值/dB');
title('双线性变换法')
五：实验结果和分析
频率/Hz
幅值
频率/Hz
相
位
归一化频率
幅值/d B
频率/hz
幅值/d B
归一化频率
幅值/d B。