哈工大威海数字信号处理实验报告.docx
《哈工大威海数字信号处理实验报告.docx》由会员分享，可在线阅读，更多相关《哈工大威海数字信号处理实验报告.docx（28页珍藏版）》请在冰豆网上搜索。
哈工大威海数字信号处理实验报告
数字信号处理
实验报告
实验日期：
2013.11.17
姓名：
何祥波
学号：
110220202
哈尔滨工业大学（威海）
实验一离散傅里叶变换的性质
一、实验目的
1、掌握离散傅里叶变换的性质，包括线性特性、时移特性、频移特性、对称性和循环卷积等性质；
2、通过编程验证傅里叶变换的性质，加强对傅里叶变换性质的认识。
二、实验原理和方法
1.线性特性
2.时移特性
3.频移特性
4.对称性
设由x（n）延拓成的周期序列为
则
共轭对称序列
共轭反对称序列
将
和
截取主周期，分别得
则
x（n）序列的实部和虚部的离散立叶变换
当x（n）为实数序列
5.循环卷积
有限长序列线性卷积与循环卷积的关系
x1（n）和x2（n）的线性卷积：
将x1（n）和x2（n）延拓成以N为周期的周期序列
则它们的周期卷积为
x1（n）和x2（n）周期延拓后的周期卷积等于他们的线性卷积的周期延拓。
三、实验内容和步骤
x1[n]，x2[n]
为长度N=8的实序列，x1[n]=[13536839]，x2[n]=[24367902]，
采用MATLAB编程验证傅里叶变换的如下性质
1.线性特性
a.给出序列x1[n]的傅里叶变换X1[k]，并画出其幅度谱和相位谱
b.给出序列x2[n]的傅里叶变换X2[k]，并画出其幅度谱和相位谱
c.给出序列Z=2*X1[k]+6*X2[k],并与序列2*x3[n]+6*x4[n]的傅里叶变换比较，
提示：
可以用fft函数实现DFT，用abs实现幅度谱，angle获得相位谱，画图用stem
程序：
（a）
clc;
clearall;
closeall;
x1=[13536839];
n=0:
length（x1）-1;
X1=fft（x1）;
figure
（1）,
subplot（2,1,1）,stem（n,abs（X1）,'r*'）;title（'X1[k]的幅度谱'）;
subplot（2,1,2）,stem（n,atan2（imag（X1）,real（X1））,'r*'）;title（'X1[k]的相位谱'）;
程序（b）
clc;
clearall;
closeall;
x2=[24367902];
n=0:
length（x2）-1;
X2=fft（x2）;
figure
（1）,
subplot（2,1,1）,stem（n,abs（X2）,'r*'）;title（'X2[k]的幅度谱'）;
subplot（2,1,2）,stem（n,atan2（imag（X2）,real（X2））,'r*'）;title（'X2[k]的相位谱'）;
程序（c）
clc;
clearall;
closeall;
x1=[13536839];
x2=[24367902];
n=0:
length（x1）-1;
X1=fft（x1）;
n=0:
length（x2）-1;
X2=fft（x2）;
Z1=2*X1+6*X2;
figure
（1）,
subplot（2,2,1）,stem（n,abs（Z1）,'r*'）;title（'Z1的幅度谱'）;
subplot（2,2,2）,stem（n,atan2（imag（Z1）,real（Z1））,'r*'）;title（'Z1的相位谱'）;
z2=2*x1+6*x2;
Z2=fft（z2）;
subplot（2,2,3）,stem（n,abs（X2）,'r*'）;title（'Z2的幅度谱'）;
subplot（2,2,4）,stem（n,atan2（imag（Z2）,real（Z2））,'r*'）;title（'Z2的相位谱'）;
2.时移特性
a.给出序列x1[n]右移3位（循环移位）后的傅里叶变换的幅度谱和相位谱，并和原始序列的幅度谱和相位谱相比较
提示：
可以用b=circshift（x1,[0,3]）实现x1序列的右移3位循环移位，或者用如下语句实现循环移位
a=1:
10%实现a=12345678910
n=3;
b1=[a（n+1:
end）,a（1:
n）]%循环移位，左移3位；实现b1=45678910123
n0=length（a）-3
b2=[a（n0+1:
end）,a（1:
n0）]%循环移位，右移3位；实现b2=89101234567
程序：
clc;
clearall;
closeall;
x1=[13536839];
b=circshift（x1,[0,3]）;
n=0:
length（x1）-1;
X1=fft（x1）;
B=fft（b）;
figure
（1）,
subplot（2,3,1）,stem（n,x1,'r*'）;title（'原始序列'）;
subplot（2,3,2）,stem（n,abs（X1）,'r*'）;title（'原始序列的幅度谱'）;
subplot（2,3,3）,stem（n,atan2（imag（X1）,real（X1））,'r*'）;title（'原始序列的相位谱'）;
subplot（2,3,4）,stem（n,b,'r*'）;title（'右移后的序列'）;
subplot（2,3,5）,stem（n,abs（B）,'r*'）;title（'右移后的序列的幅度谱'）;
subplot（2,3,6）,stem（n,atan2（imag（B）,real（B））,'r*'）;title（'右移后的序列的相位谱'）;
3.对称性
（1）利用x1[n]构造圆周共轭对称序列和圆周共轭反对称序列（第三章学习），讨论如下问题
（a）画出圆周共轭对称序列的傅里叶变换的幅度谱和相位谱
（b）画出圆周共轭对称序列的傅里叶变换的实部和虚部
（c）利用x1[n]构造共轭对称序列和共轭反对称序列（第二章学习），讨论与圆周共轭对称圆周共轭反对称的区别
提示：
应先构造出圆周共轭对称或共轭反对称序列。
可以用
x1=[13536839]
x2=fliplr（x1）%翻褶序列，x2=[93863531]
x3=circshift（x2,[1,1]）%x2的右移一位序列，注意是循环移位，x3=[19386353]
xe=（x1+conj（x3））/2%构造出圆周共轭对称序列
xo=（x1-conj（x3））/2%构造出圆周共轭反对称序列
函数real、imag、abs、angle分别为求实部、虚部、幅度和相位。
思考：
共轭对称和共轭反对称序列如何构造？
程序（a）：
clc;
clearall;
closeall;
x1=[13536839]
x2=fliplr（x1）%翻褶序列，x2=[93863531]
x3=circshift（x2,[1,1]）%x2的右移一位序列，注意是循环移位，x3=[19386353]
xe=（x1+conj（x3））/2%构造出圆周共轭对称序列
n=0:
length（xe）-1;
Xe=fft（xe）;
figure
（1）,
subplot（2,1,1）,stem（n,abs（Xe）,'r*'）;title（'圆周共轭对称的傅里叶变换的幅度谱'）;
subplot（2,1,2）,stem（n,atan2（imag（Xe）,real（Xe））,'r*'）;title（'圆周共轭对称的傅里叶变换的相位谱'）;
程序（b）：
clc;
clearall;
closeall;
x1=[13536839]
x2=fliplr（x1）%翻褶序列，x2=[93863531]
x3=circshift（x2,[1,1]）%x2的右移一位序列，注意是循环移位，x3=[19386353]
xe=（x1+conj（x3））/2%构造出圆周共轭对称序列
n=0:
length（xe）-1;
Xe=fft（xe）;
figure
（1）,
subplot（2,1,1）,stem（n,real（Xe）,'r*'）;title（'圆周共轭对称序列的傅里叶变换的实部'）;
subplot（2,1,2）,stem（n,imag（Xe）,'r*'）;title（'圆周共轭对称序列的傅里叶变换的虚部'）;
程序（c）：
clc;
clearall;
closeall;
x1=[13536839];
N=length（x1）;
n_axis=[0:
N-1];
x1N=zeros（1,N）;
x1N（0+1）=x1（0+1）;
fornum=0:
N-1
ifN-num==N
index=1;
else
index=（N-num）+1;
end
x1N（num+1）=x1（index）;
end
x1_e=1/2*（x1+x1N）;
x1_o=1/2*（x1-x1N）;
figure
（1）;
subplot（3,1,1）;stem（n_axis,x1）;title（['原序列']）;
subplot（3,1,2）;stem（n_axis,x1_e）;title（'共轭对称部分'）;
subplot（3,1,3）;
stem（n_axis,x1_o）;title（'共轭反对称部分'）;
4.循环卷积与线性卷积
a.计算序列x1[n]和x2[n]的线性卷积y[n]
b.计算x1[n]和x2[n]的傅里叶变换X1[k]和X2[k]，Y[k]=X1[k]*X2[k]，求Y[k]的反傅里叶变换y2[n]（y2[n]即为x1[n]和x2[n]的循环卷积），比较y[n]与y2[n]的异同.
程序（a）：
clc;
clearall;
closeall;
x1=[13536839];
x2=[24367902];
y=conv（x1,x2）;
n=0:
length（y）-1;
subplot（1,1,1）;stem（n,y）;title（'线性卷积'）;axis（[0,17,0,200]）;
程序（b）:
clc;
clearall;
closeall;
x1=[13536839];
x2=[24367902];
y=conv（x1,x2）;
X1=fft（x1）;
X2=fft（x2）;
Y2=X1.*X2;
y2=ifft（Y2）;
n1=0:
length（x1）-1;
n2=0:
length（x2）-1;
n3=0:
length（y）-1;
n4=0:
length（y2）-1;
subplot（2,2,1）;stem（n1,x1）;title（'序列x1[n]'）;
subplot（2,2,2）;stem（n2,x2）;title（'序列x2[n]'）;
subplot（2,2,3）;stem（n3,y）;title（'序列x1[n]和x2[n]的循环卷积y[n]'）;
subplot（2,2,4）;stem（n4,y2）;title（'Y[k]的反傅里叶变换y2[n]'）;
5.补零
用MATLAB计算如下N点序列的M点DFT：
（1）取N=8，M=8
（2）取N=8，M=16
（3）取N=8，M=32
根据实验结果，分析同一序列不同点数的离散傅里叶变换的特点。
程序：
clc;
clearall;
closeall;
N=8;
M=[81632];
x=ones（1,N）;
k_axis=[0:
max（M）-1];
fft_x8=fft（x,M
（1））;fft_x8=fftshift（fft_x8）;
fft_x16=fft（x,M
（2））;fft_x16=fftshift（fft_x16）;
fft_x32=fft（x,M（3））;fft_x32=fftshift（fft_x32）;
figure
（1）;
subplot（3,1,1）;stem（k_axis（1:
M
（1））,abs（fft_x8））;title（'8点序列的8点DFT'）;
subplot（3,1,2）;stem（k_axis（1:
M
（2））,abs（fft_x16））;title（'8点序列的16点DFT'）;
subplot（3,1,3）;stem（k_axis（1:
M（3））,abs（fft_x32））;title（'8点序列的32点DFT'）
实验二时域采样与频域采样
1.时域采样定理的验证，分析采样序列的特性。
对xa（t）=Ae-atsin（Ω0t）u（t）进行采样，可得到采样序列
xa（n）=xa（nT）=Ae-anTsin（Ω0nT）u（n）,
其中A为幅度因子，a为衰减因子，Ω0是模拟角频率，T为采样间隔。
a.取采样频率fs=1kHz,即T=1ms，观察|X（ejω）|。
b.改变采样频率，fs=300Hz，观察|X（ejω）|的变化；
c.进一步降低采样频率，fs=200Hz，观察频谱混叠是否明显存在，说明原因，并记录|X（ejω）|曲线。
2.频域采样定理的验证。
给定信号
编写程序分别对频谱函数
在区间
上等间隔采样32点和16点，得到X32（k）和X16（k），再分别对其进行32点和16点的IFFT，得到x32（n）和x16（n），分别画出
、X32（k）和X16（k）的幅度谱。
并绘图显示x（n）、x32（n）和x16（n）的图形，进行对比分析，验证频域采样定理。
程序
（1）：
%采样
A=44.128;
a=50*sqrt
（2）.*pi
W0=50*sqrt
（2）.*pi
Fs1=1000;%采样频率Hz
Fs2=300;
Fs3=200;
Ts1=1/Fs1;
Ts2=1/Fs2;
Ts3=1/Fs3;
T1=0.1;%采样时间为0.1s
figure
（1）;
t1=[0:
Ts1:
T1-Ts1]%采样频率1000
x1=A.*exp（-a*t1）.*sin（W0*t1）;
subplot（3,1,1）;
stem（t1,x1,'.'）;
t2=[0:
Ts2:
T1-Ts2]%采样频率300
x2=A.*exp（-a*t2）.*sin（W0*t2）;
subplot（3,1,2）;
stem（t2,x2,'.'）;
t3=[0:
Ts3:
T1-Ts3]%采样频率200
x3=A.*exp（-a*t3）.*sin（W0*t3）;
subplot（3,1,3）;
stem（t3,x3,'.'）;
figure
（2）;
X1=fft（x1）;
AMP1=abs（X1）;
PHA1=angle（X1）;
subplot（3,1,1）
plot（t1,AMP1）
gridon
X2=fft（x2）;
AMP2=abs（X2）;
PHA2=angle（X2）;
subplot（3,1,2）
plot（t2,AMP2）
gridon
X3=fft（x3）;
AMP3=abs（X3）;
PHA3=angle（X3）;
subplot（3,1,3）
plot（t3,AMP3）
gridon
程序
（2）：
M=27;
N=32;
n=0:
M;
xa=0:
floor（M/2）;
xb=ceil（M/2）-1:
-1:
0;
xn=[xa,xb];
Xk=fft（xn,1024）;
X32k=fft（xn,32）;
x32n=ifft（X32k）;
X16k=X32k（1:
2:
N）
x16n=ifft（X16k,N/2）;
subplot（3,2,2）;
stem（n,xn,'.'）;
boxon
title（'（b）三角波序列x（n）'）
xlabel（'n'）;
ylabel（'x（n）'）;axis（[0,32,2,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）'.'）;
boxon
title（'（c）16点频域采样'）
xlabel（'k'）;
ylabel（'|X_1_6（k）|'）;axis（[0,8,0,20]）;
n1=0:
N/2-1;
subplot（3,2,4）;
stem（n1,x16n,'.'）;
boxon
title（'（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）,'.'）;
boxon
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,'.'）;
boxon
title（'（f）32点IDFT[X_3_2（k）]'）;
xlabel（'n'）;
ylabel（'x_3_2（n）'）;axis（[0,32,0,20]）;
实验三IIR滤波器设计
一、实验目的
1、掌握冲激响应不变法和双线性变换法设计IIR滤波器的原理及具体设计方法，熟悉用双线性设计法设计低通IIR数字滤波器的计算机程序；
2、熟悉模拟Butterworth滤波器的设计，掌握冲激响应不变法和双线性变换法设计数字IIR滤波器的方法。
二、实验原理和方法
IIR滤波器设计的过程可以首先设计模拟滤波器，然后采用冲激响应不变法和双线性变换法设计IIR数字滤波器。
Butterworth滤波器
，其中
为3dB截止频率，N为滤波器阶次，均待定。
模拟滤波器的设计步骤：
首先根据数字滤波器设计要求计算模拟滤波器指标；其次要求出滤波器的阶次N和
三、实验内容
1.采样频率为1Hz，设计一个Butterworth低通数字滤波器，其中通带临界频率
，通带内衰减小于1dB（
），阻带临界频率
，阻带内衰减大于25dB（
）。
求这个数字滤波器的传递函数H（z），输出它的幅频特性曲线。
利用冲激响应不变法和双线性变换法实现该滤波器，并将结果进行比较。
程序：
clc;
clearall;
closeall;
J=sqrt（-1）;
freq_p=0.2*pi;%通带截止频率0.2pi
alpha_p=1;%通带最大衰减1dB
freq_s=0.3*pi;%阻带截止频率0.3pi
alpha_s=25;%阻带最小衰减15dB
Fs_array=[12];%采样频率
Color_array=['r','k'];
forIII=1:
2
color_used=Color_array（III）;
Fs=Fs_array（III）;
T=1/Fs;
analog_freq_p=freq_p/T;%数字频率w到模拟角频率W的转换W=w/T;
analog_freq_s=freq_s/T;
%根据模拟滤波器的指标求巴特沃兹滤波器的阶次
k_sp=sqrt（10^（alpha_p/10）-1）/sqrt（10^（alpha_s/10）-1）;%求ksp
lamda_sp=2*pi*analog_freq_s/（2*pi*analog_freq_p）;%求lanmdasp
N=ceil（-log10（k_sp）/log10（lamda_sp））;%求N
Wc=analog_freq_p*（10^（analog_freq_p/10）-1）^（-1/2/N）;%3dB截至频率
p=zeros（1,N）;
forI=1:
N
p（I）=exp（J*pi/2+J*pi*（2*I-1）/2/N）;%求传输函数极点
end
%根据得到的零点、极点得到归一化的传输函数H（p）
z=[];%表示无零点
[ba]=zp2tf（z,p,1）;%得到传递函数tf:
trans-function,第三个参数为增益=1
%根据p=s/Wc把归一化传输函数H（p）的系数转换为实际的模拟传输函数H（s）
bs=zeros（1,length（b））;%H（s）的分子系数
as=zeros（1,length（a））;%对与N阶模拟滤波器，设其传输函数分子分母系数分别为bs,as
temp=[N:
-1:
0];%以下为对H（s）进行去归一化
temp=Wc.^temp;
bs=b./temp;
as=a./temp;
temp=as
（1）;
as=as/temp;%验证与书上结果是否一致
bs=bs/temp;
freq_analog=linspace（0,2*analog_freq_s）;
H_analog=freqs（bs,as,freq_analog）;
figure
（1）;
plot（freq_analog,20*log10（abs（H_analog））,color_used）;gridon;holdon;
xlabel（'模拟角频率（rad/s）'）;ylabel（'幅频响应（dB）'）;
title（'模拟滤波器的频率响应'）;
%根据H（s）获得H（z），使用脉冲响应不变法
[bz,az]=impinvar（bs,as,Fs）;%根据模拟滤波器系数b,a，按照采样率fs转换为系数为bz,az的数字IIR滤波器
freq_digital=[0:
0.1:
2*pi-0.1];
H=freqz（bz,az,freq_digital）;
figure
（2）;
plot（freq_digital/pi,20*log10（abs（H））,color_used）;gridon;holdon;
xlabel（'数字频率/pi'）;ylabel（'幅频响应（dB）'）;
title（'数字IIR滤波器频率响应'）;
end
figure
（1）,legend（['Fs='num2str（Fs_array
（1））'Hz'],['Fs
升级会员 