实验三DFT和DCT及频域滤波.docx
《实验三DFT和DCT及频域滤波.docx》由会员分享，可在线阅读，更多相关《实验三DFT和DCT及频域滤波.docx（21页珍藏版）》请在冰豆网上搜索。
实验三DFT和DCT及频域滤波
一．
实验名称：
数字信号的DFT/DCT及频域滤波
二．实验目的
1．熟练掌握数字信号（1D）及数字图像（2D）离散傅立叶变换（DFT）及离散余弦变换（DCT）方法、基本原理及实现流程。
熟悉两种变换的性质，并能对DFT及DCT的结果进行必要解释。
2．深入理解离散信号采样频率、奈奎斯特频率及频率分辨率等基本概念，弄清它们之间的相互关系。
了解离散傅里叶变换（DFT）中频率泄露的原因，以及如何尽量减少频率泄露影响的途径。
3．熟悉和掌握利用MATLAB工具进行1D/2DFFT及DCT的基本步骤、MATLAB函数使用及对具体变换的处理流程。
4．能熟练应用MATLAB工具对数字图像进行FFT及DCT处理，并能根据需要进行必要的频谱分析和可视化显示。
三．实验原理
1、傅立叶变换
●傅立叶变换：
非周期函数表示为正弦和/或余弦乘以加权函数的积分。
●一维连续Fourier变换
对函数f（x）进行傅立叶变换得到F（u）
逆变换，即将F（u）变换到f（x）为
●一维离散Fourier变换
正变换（DFT）
逆变换（IDFT）
●用幅值和相位表示傅立叶变换
2、离散余弦变换
●1D-DCT
●IDCT变换
●矩阵形式
四．实验步骤
1.1D数字信号的FFT及频谱分析
给定如下式
（1）所示的1D连续信号：
1）设采样频率
=1000Hz，对信号
进行离散化，并画出一个周内的信号振幅随时间变化的波形图。
2）对离散信号
进行傅立叶变换，分别画出频谱中心化及有效频率范围（不含负频）2种方式下的幅值（
）随频率
变化的分布图，要求纵横坐标正确标注物理量和单位。
3）对式
（1）信号，加随机噪声，重复步骤
（1）和
（2）的处理过程。
4）通过对变换结果的分析，说明采样频率
、奈奎斯特（Nyquist）频率（
）及采样时间间隔△T三者之间的相互关系，并简要描述模拟信号的采样定理。
图1
注：
加随机噪声的处理只需改变初始信号即可，其它步骤同理。
2.数字音频信号的DFT
1）读取一段0.5s的预先录制的数字音频信号（“yes.wav”或“no.wav”文件中任选其一），画出随时间变化的声波波形图。
2）对数字音频信号进行离散傅立叶变换（DFT），分别画出频谱中心化及有效频率范围（不含负频）2种方式下的幅值（
）随频率
变化的分布图，要求纵横坐标正确标注物理量和单位。
图2
3.数字音频信号的DCT和IDCT
1）对上述音频信号做离散余弦变换（DCT），画出DCT变换系数（变换结果）图，并对变换结果进行必要的解释，说明DCT变换的主要用途。
要求按DCT原理自行编写实现代码，不允许直接调用MATLAB的dct（）函数。
2）按原理自行一段MATLAB代码，对第
（1）步处理结果进行离散余弦反变换（IDCT），将计算结果与原始音频信号进行比较，检验编写代码的正确性。
3）编写一段MATLAB代码，利用快速傅立叶变换（FFT）程序实现快速DCT算
法（FCT），并将计算结果与直接调用dct（）的处理结果进行比较，检验编写代码的正确性。
图3
4.综合应用题：
实际信号的频谱分析及频域滤波
1）编写一从保存在本地磁盘的文本文件中读入一实际数字信号，已知该信号的时间采样率为dt=2ms。
文件中的信号由301个等长的按列排列的一维列信号组成，每个一维列信号有251个采样点，信号实际计时起点为1800ms，延时长度为L=（251-1）*2ms=500ms。
请读出其中的某一列信号，并画出该信号振幅随时间变化的波形图，以ms为时间单位。
2）对第一步中抽取的其中一列信号做快速傅里叶变换（FFT），分别画出频谱中心化的对称频谱和只含有正半抽的信号频谱图，并对该信号做简要的频谱分析。
要求规范的标注纵横坐标实际物理量和对应的单位。
3）设定截止频率D0=100，试在同一张图上以不同线型画出n=1，2，4阶下的巴特沃思（Butterworth）低通滤波器（一维）的频率响应曲线。
要求标注规范地纵横坐标实际物理量和对应的单位。
4）选择合适的D0，利用上述2阶Butterworth低通滤波器，对第
（1）步读取的列信号进行滤波实验。
并分析截止频率对滤波效果的影响。
图4
五．实验结果及分析
1.
图5
1D数字信号的FFT及频谱分析
●
图6
分析：
通过对变换结果的分析采样频率
、奈奎斯特频率
及采样时间间隔三者之间有如下的相互关系
●模拟信号的采样定理：
当采样频率
大于连续信号最高频率
时，采样之后的数字信号完整地保留了原始信号中的信息，是连续信号离散化的基本依据。
2.
图7
数字音频信号的DFT
●分析：
从上面的频谱可以看出该音频的频率成分主要集中在700Hz以内，共有4支峰。
3.
图8
数字音频信号的DCT和IDCT
●分析：
DCT变换将信号从时域转换到变换域上，通过对变换后的系数分析，原能量集中在少数系数上，可以提高编码效率，压缩数据。
4.
图9
综合应用题：
实际信号的频谱分析及频域滤波
图10
图11
图12
●该信号的第二列数据的有效频率主要集中在0-100Hz以内
●通过设置不同截止频率的Butterworth低通滤波器，可以看出对于低频信号，截止频率从100Hz减小，原信号的有效频率则被滤去的变多，若是图像信号，则表现出图像变的模糊。
六.实验心得体会和建议
●心得体会：
通过这次实验使我深刻了解了奈奎斯特定理、DFT和DCT的基本原理以及巴特沃斯低通滤波器的构造，与此同时在上机实验中熟悉了MATLAB编写FFT和DCT的基本方法及步骤。
在对相关频谱的分析过程中更加深刻的理解了原理，及相关用途。
●建议：
可以让大家不调用FFT函数，直接编写DCT。
七.程序源代码
1.1D数字信号的FFT及频谱分析
fs=1000;
N=256;
n=0:
N-1;
t=n/fs;
x=2*sin（30*pi*t）+0.5*cos（120*pi*t）+4*sin（240*pi*t）;%输入信号
plot（t,x）,xlim（[0,1/15]）;%画出一个周期内的信号振幅随时间变化的波形图
%====================================================================
y=fft（x）;
mag=abs（y）;
y=fftshift（y）;%频谱中心化
mag0=abs（y）;
%====================================================================
M=length（y）;%fft频率轴点数（maybedifferent）
f=（0:
M-1）*fs/M;%频率采样序列（矢量）
%====================================================================
fchar=num2str（fs）;%采样率转化为char（本文）
nchar=num2str（N）;%样点数转化为char（本文）
ltext=strcat（'fs=',fchar,'Hz',',N=',nchar,'points'）;%拼title字符
%====================================================================
subplot（311）,plot（t,x）;%随时间变化的振幅
xlabel（'t/s'）;ylabel（'振幅'）;xlim（[0,0.1]）;
title（['x=2*sin（30*pi*t）+0.5*cos（120*pi*t）+4*sin（240*pi*t）']）;
gridon;
f0=f-f（M/2）;
subplot（312）,plot（f0,（mag0）*2/N）;%随频率变化的振幅
xlabel（'频率/Hz'）;ylabel（'振幅'）;
title（['全部频率：
',ltext]）;gridon;
subplot（313）,plot（f（1:
M/2）,（mag（1:
M/2））*2/N）;%绘制有效频谱
xlabel（'频率/Hz'）;ylabel（'振幅'）;
title（['有效频率：
',ltext]）;gridon;
注：
对于加性噪声只需将x变为
x=2*sin（30*pi*t）+0.5*cos（120*pi*t）+4*sin（240*pi*t）+randn（size（t））;
即可，其它程序同理。
2.数字音频信号的DFT
[x,fs]=wavread（'yes.wav'）;
N=4000;
n=0:
N-1;
t=n/fs;
plot（t,x）,%画出一个周期内的信号振幅随时间变化的波形图
%====================================================================
y=fft（x）;
mag=abs（y）;
y=fftshift（y）;%频谱中心化
mag0=abs（y）;
%====================================================================
M=length（y）;%fft频率轴点数（maybedifferent）
f=（0:
M-1）*fs/M;%频率采样序列（矢量）
%====================================================================
fchar=num2str（fs）;%采样率转化为char（本文）
nchar=num2str（N）;%样点数转化为char（本文）
ltext=strcat（'fs=',fchar,'Hz',',N=',nchar,'points'）;%拼title字符
%====================================================================
subplot（311）,plot（t,x）;%随时间变化的振幅
xlabel（'t/s'）;ylabel（'振幅'）;xlim（[0,0.5]）;
title（['yes']）;gridon;
f0=f-f（M/2）;
subplot（312）,plot（f0,（mag0）*2/N）;%随频率变化的振幅
xlabel（'频率/Hz'）;ylabel（'振幅'）;
title（['全部频率：
',ltext]）;gridon;
subplot（313）,plot（f（1:
M/2）,（mag（1:
M/2））*2/N）;%绘制有效频谱
xlabel（'频率/Hz'）;ylabel（'振幅'）;
title（['有效频率：
',ltext]）;gridon;
3.数字音频信号的DCT和IDCT
[y,Fs]=wavread（'yes.wav'）;
N=length（y）;
t=（0:
N-1）/Fs;
subplot（311）,plot（t,y）;
xlabel（'t/s'）;ylabel（'振幅'）;
title（'yes波形'）;
G=zeros（N,N）;
forx=1:
N;
foru=2:
N;
G（x,u）=cos（（2*（x-1）+1）*（u-1）*pi/（2*N））;
end
end
G（1,:
）=sqrt（2/N）;
G（1,1）=sqrt（1/N）;
F=G*y;
subplot（312）,plot（（1:
N）/Fs,F）
title（'DCTcoeffcients'）;
xlabel（'t/s'）;ylabel（'振幅'）;
f=G'*F;
subplot（313）,plot（t,f/Fs*4）;
title（'TherecoveryofthesignalfromDCTcoefficients'）
xlabel（'t/s'）;ylabel（'振幅'）;
4.综合应用题：
实际信号的频谱分析及频域滤波
clc,clear,closeall;
matrix=importdata（'seismic_nsamp251_tr301_2ms.txt'）;
[p,q]=size（matrix）;
fs=500;
r=2;%这里取第2列
dt=1/fs;
L=（p-1）*dt;%每列的延时长度
t0=1.8;
to_r=t0+L*（r-1）;%第r列的开始时间
sl=matrix（:
r）;%取出第r列的数据
n=0:
p-1;
t=to_r+n/fs;
%t=n/fs;
%====================================================================
%划出振幅随时间的变化曲线
%====================================================================
figure
（1）;
plot（t,sl'）;%随时间变化的振幅
xlabel（'t/ms'）;ylabel（'振幅'）;
title（[strcat（'第',num2str（r）,'列'）]）;gridon;
%====================================================================
y=fft（sl,p-1）;
mag=abs（y）;
y=fftshift（y）;%频谱中心化
mag0=abs（y）;
%====================================================================
M=length（y）;%fft频率轴点数（maybedifferent）
f=（0:
M-1）*fs/M;%频率采样序列（矢量）
%====================================================================
%画出中心频率图
f0=f-f（M/2）;
figure;
subplot（211）,plot（f0,mag0）;%随频率变化的振幅
xlabel（'频率/Hz'）;ylabel（'振幅'）;
title（['全部频率']）;gridon;
%====================================================================
subplot（212）,plot（f（1:
M/2）,mag（1:
M/2））;%绘制有效频谱
xlabel（'频率/Hz'）;ylabel（'振幅'）;
title（['有效频率']）;gridon;
%===============================================================
%生成Butterworth矩阵
H=zeros（3,M）;
JS=[124];
D0=100;
forn=1:
M
H（1,n）=1/（1+（（n-（M）/2）/D0）^（2*JS
（1）））;%一阶Butterworth矩阵
H（2,n）=1/（1+（（n-（M）/2）/D0）^（2*JS
（2）））;%二阶Butterworth矩阵
H（3,n）=1/（1+（（n-（M）/2）/D0）^（2*JS（3）））;%四阶Butterworth矩阵
end
%===============================================================
%在同一图中画出butterworth函数
figure（3）;
subplot（311）,plot（f0,abs（H（1,:
））,'-b'）;holdon;
subplot（312）,plot（f0,abs（H（2,:
））,'-b'）;holdon;
subplot（313）,plot（f0,abs（H（3,:
））,'-b'）;holdon;
xlim（[f0
（1）,f0（M）]）;xlabel（'f0/Hz'）;
%=============================================================
z1=y'.*H（2,:
）;
figure（4）,plot（f0,abs（z1））,gridon;
八.思考题
1.分别阐述和解释什么叫信号的采样频率、奈奎斯特（Nyquist）频率、时间采样率及频率分辨率？
a）信号的采样频率：
对模拟信号进行A/D采样时，每秒钟对信号采样的点数。
b）奈奎斯特频率：
指最低允许的抽样率，是带限信号频率宽度的2倍。
c）采样时间：
每进行一次采样所需要的时间。
d）频率分辨率：
指将两个相邻谱峰分开的能力，也指分辨两个不同频率信号的最小间隔。
2.根据所学知识，简要叙述离散傅立叶变换（DFT）和离散余弦变换（DCT）在数字信号处理中的主要用途。
主要用途：
a）DFT：
频域滤波，图像增强，频谱分析。
b）DCT：
数据压缩，提高编码效率。
3.简要叙述频率滤波与时域滤波在处理上有什么不同。
区别：
对于不规则的噪声值，通过FFT进行频率滤波可以得到很好的平滑效果，但是FFT不可避免的会存在频谱泄露；而时域在边缘提取上，要比频域的处理更优秀。
4.实验内容1中，式
（1）表示的信号是平稳信号还是非平稳信号，简要说明平稳信号和非平稳信号的区别。
a）式
（1）表示的信号是平稳信号
b）区别：
通过对统计特征的对比，非平稳信号的统计特性（均值、方差等）随着时间变化而变化，而平稳信号的统计特性不随时间变化。
报告评分：
指导教师签字：
升级会员 