数字信号的 DFTDCT 及频域滤波.docx
《数字信号的 DFTDCT 及频域滤波.docx》由会员分享，可在线阅读，更多相关《数字信号的 DFTDCT 及频域滤波.docx（16页珍藏版）》请在冰豆网上搜索。
数字信号的DFTDCT及频域滤波
实验报告
实验课程:
光电图像处理
姓名：
学号：
实验地点:
信软楼309
指导老师：
实验时间：
2016年3月30日
一．实验名称:
数字信号的DFT/DCT及频域滤波
二．实验目的
1、熟练掌握数字信号（1D）及数字图像（2D）离散傅立叶变换（DFT）及离散余弦变换（DCT）方法、基本原理及实现流程。
熟悉两种变换的性质，并能对DFT及DCT的结果进行必要解释。
2、深入理解离散信号采样频率、奈奎斯特频率及频率分辨率等基本概念，弄清它们之间的相互关系。
了解离散傅里叶变换（DFT）中频率泄露的原因，以及何尽量减少频率泄露影响的途径。
3、熟悉和掌握利用MATLAB工具进行1D/2DFFT及DCT的基本步骤、MATLAB函数使用及对具体变换的处理流程。
4、能熟练应用MATLAB工具对数字图像进行FFT及DCT处理，并能根据需要进行必要的频谱分析和可视化显示。
三.实验原理
一、一维连续Fourier变换
二、离散余弦变换（DCT）
DFT的特殊形式，被展开的函数是实偶函数的傅立叶变换，即只有余弦项。
变换核固定，利于硬件实现。
可分离性，一次二维DCT可分解为两次一维DCT。
其有快速算法FCT（可由FFT计算得到）。
1D-DCT
矩阵形式：
F=Gf，f=GTF
2D-IDCT
三、巴特沃思（Butterworth）低通滤波器
BLPF传递函数
这里D0为截止频率，n为阶数。
当D（u,v）=D0时,H（u,v）降为最大值的50%。
P、Q意义如前所述
四．实验步骤
1、1D数字信号的FFT及频谱分析
给定如下式
（1）所示的1D连续信号：
（1）设采样频率fs=1000Hz，对信号x进行离散化，并画出一个周期内的信号振幅随时间变化的波形图。
（2）对离散信号x进行傅立叶变换，分别画出频谱中心化及有效频率范围（不含负频）2种方式下的幅值（|F（u）|）随频率（u）变化的分布图，要求纵横坐标正确标注物理量和单位。
（3）对式
（1）信号，加随机噪声，重复步骤
（1）和
（2）的处理过程。
（4）通过对变换结果的分析，说明采样频率fs、奈奎斯特（Nyquist）频率（fn）及采样时间间隔△T三者之间的相互关系，并简要描述模拟信号的采样定理
2、数字音频信号的DFT
（1）读取一段0.5s的预先录制的数字音频信号（“yes.wav”或“no.wav”文件中任选其一），画出随时间变化的声波波形图，如图3-1所示。
（2）对数字音频信号进行离散傅立叶变换（DFT），分别画出频谱中心化及有效频率范围（不含负频）2种方式下的幅值（|F（u）|）随频率（u）变化的分布图，要求纵横坐标正确标注物理量和单位。
3、数字音频信号的DCT和IDCT
（1）对图3-1中的音频信号做离散余弦变换（DCT），画出DCT变换系数（变换结果）图，并对变换结果进行必要的解释，说明DCT变换的主要用途。
要求按DCT原理自行编写实现代码，不允许直接调用MATLAB的dct（）函数。
（2）按原理自行一段MATLAB代码，对第
（1）步处理结果进行离散余弦反变换（IDCT），将计算结果与原始音频信号进行比较，检验编写代码的正确性。
（3）编写一段MATLAB代码，利用快速傅立叶变换（FFT）程序实现快速DCT算
法（FCT），并将计算结果与直接调用dct（）的处理结果进行比较，检验编写代码的正确性。
（已布置课程作业，这里作为选作实验！
）
4、综合应用题：
实际信号的频谱分析及频域滤波
（1）从保存在本地磁盘的文本文件中读入一实际数字信号，该磁盘文件名为:
“seismic.txt”，已知该信号的时间采样率为dt=2ms。
文件中的信号由301
个等长的按列排列的一维列信号组成，每个一维列信号有251个采样点，信号实际计时起点为1800ms，延时长度为L=（251-1）*2ms=500ms。
请读出其中的某一列信号，并画出该信号振幅随时间变化的波形图，以ms为时间单位。
（2）对第一步中抽取的其中一列信号做快速傅里叶变换（FFT），分别画出频谱中心化的对称频谱和只含有正半抽的信号频谱图，并对该信号做简要的频谱分析。
要求规范的标注纵横坐标实际物理量和对应的单位。
（3）设定截止频率D0=100，试在同一张图上以不同线型画出n=1，2，4阶下的巴特沃思（Butterworth）低通滤波器（一维）的频率响应曲线。
要求标注规范地纵横坐标实际物理量和对应的单位。
其中，D0为截止频率，N为滤波器长度，n为滤波器阶数。
（4）选择合适的D0，利用上述2阶Butterworth低通滤波器，对第
（1）步读取的列信号进行滤波实验。
并分析截止频率对滤波效果的影响。
五．实验结果及分析
1、1D数字信号的FFT及频谱分析
加噪声之后：
结果分析：
对比两张图可以看多唯一的区别是加噪声图片在频率为0Hz左右振幅变大。
二、数字音频信号的DFT
结果如下：
3、数字音频信号的DCT和IDCT
结果分析：
从图中可以看出编写的一维DCT函数与直接调用dct（）函数效果几乎一样。
4、综合应用题：
实际信号的频谱分析及频域滤波
Butterworth低通滤波器（一维）的频率响应曲线
2阶Butterworth低通滤波器对数据进行分析：
从图中可以看出截止频率小的滤波器滤波效果要好。
六.实验心得体会和建议
（1）、直接调用滤波函数的效果不是很好。
很多好的函数还是得自己写出来才行。
（2）、实验开始的时候要看清楚给的数据需要用的是一味CDT还是二维DCT变换。
因为我在看是的时候用的是二维的DCT变换结果得到的结果差不多。
这也是我的疑问。
七.程序源代码
1、1D数字信号的FFT及频谱分析
clc;clf;clearall;closeall;
fs=1000;%采样频率（设定,Hz）
N=128;%采样点数（可任意设定）,Suchas64,128,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）+rand（）;
%多分量信号
subplot（411）,plot（t,x,'.'）;%原始信号
title（['原始信号:
','N=',num2str（N）,'points']）;
xlabel（'时间/s'）;ylabel（'振幅'）;
xlim（[0max（t）]）;
y=fft（x,512）;%对输入信号的FFT
mag=abs（y）;%FFT幅值（未中心化）
y=fftshift（y）;%频谱中心化,matlab内部函数
mag0=abs（y）;%FFT幅度（中心化）
M=length（y）;%FFT频率轴点数（maybedifferent）
df=fs/（M-1）;%频率步长
f=（0:
M-1）*df;%频率采样序列（矢量）
fchar=num2str（fs）;%采样率转化为char（本文）
nchar=num2str（M）;%样点数转化为char（本文）
ltext=strcat（'fs=',fchar,'Hz',',N=',nchar,'points'）;%拼接title字符
subplot（412）,plot（f,mag）;%随频率变化的振幅
xlabel（'频率/Hz'）;ylabel（'振幅'）;
title（['全部频率：
',ltext]）;gridon;
f0=f-f（M/2+1）;
subplot（413）,plot（f0,mag0）;%随频率变化的振幅
xlabel（'频率/Hz'）;ylabel（'振幅'）;
title（['对称频率：
',ltext]）;gridon;
xlim（[-fs/2,fs/2]）;
subplot（414）,plot（f（1:
M/2）,mag（1:
M/2））;%绘制有效频谱
xlabel（'频率/Hz'）;ylabel（'振幅'）;
title（['有效频率：
',ltext]）;gridon;
xlim（[0,fs/2]）;
2、数字音频信号的DFT
clc,clf,clearall,closeall;
fs=1000;%采样频率（设定,Hz）
[x,fs]=audioread（'yes.wav'）;
N=length（x）;%采样点数（可任意设定）,Suchas64,128,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）;%多分量信号
subplot（411）,plot（t,x）;%原始信号
title（['原始信号:
','N=',num2str（N）,'points']）;
xlabel（'时间/s'）;ylabel（'振幅'）;
xlim（[0max（t）]）;
y=fft（x,N）;%对输入信号的FFT
mag=abs（y）;%FFT幅值（未中心化）
y=fftshift（y）;%频谱中心化,matlab内部函数
mag0=abs（y）;%FFT幅度（中心化）
M=length（y）;%FFT频率轴点数（maybedifferent）
df=fs/（M-1）;%频率步长
f=（0:
M-1）*df;%频率采样序列（矢量）
fchar=num2str（fs）;%采样率转化为char（本文）
nchar=num2str（M）;%样点数转化为char（本文）
ltext=strcat（'fs=',fchar,'Hz',',N=',nchar,'points'）;%拼接title字符
subplot（412）,plot（f,mag）;%随频率变化的振幅
xlabel（'频率/Hz'）;ylabel（'振幅'）;
title（['全部频率：
',ltext]）;gridon;
f0=f-f（M/2+1）;
subplot（413）,plot（f0,mag0）;%随频率变化的振幅
xlabel（'频率/Hz'）;ylabel（'振幅'）;
title（['对称频率：
',ltext]）;gridon;
xlim（[-fs/2,fs/2]）;
subplot（414）,plot（f（1:
M/2）,mag（1:
M/2））;%绘制有效频谱
xlabel（'频率/Hz'）;ylabel（'振幅'）;
title（['有效频率：
',ltext]）;gridon;
xlim（[0,fs/2]）;
3.数字音频信号的DCT和IDCT
clc,clf,clear,close;
fs=1000;%采样频率（设定,Hz）
[X,fs]=audioread（'yes.wav'）;
[m,n]=size（X）;
AM=zeros（m,m）;
AN=zeros（n,n）;
fori=0:
m-1
forj=0:
m-1
ifi==0
AM（i+1,j+1）=sqrt（1/m）*cos（（（2*j+1）*i*pi）/（2*m））;
else
AM（i+1,j+1）=sqrt（2/m）*cos（（（2*j+1）*i*pi）/（2*m））;
end
end
end
X=double（X）;%注意变换后的矩阵数据类型为double
Y=AM*X;
Y1=AM'*Y;%实现DCT逆变换
M=length（Y）;%FFT频率轴点数
df=fs/（M-1）;%频率步长
f=（0:
M-1）*df;
subplot（3,1,1）
plot（f,Y）;
xlabel（'频率/Hz'）;ylabel（'振幅'）;
title（'DCT图像'）;
subplot（3,1,2）
plot（f,Y1）;
xlabel（'频率/Hz'）;ylabel（'振幅'）;
title（'IDCT的图像'）;
Y2=dct（X）;
subplot（3,1,3）
plot（f,Y2）;
xlabel（'频率/Hz'）;ylabel（'振幅'）;
title（'直接调用dct（）函数的图像'）;
4、综合应用题：
实际信号的频谱分析及频域滤波
clc,clf,clear,close;
I=importdata（'seismic.txt'）;
y=I（:
2）;
dt=2;
fs=1/0.002;
L=251;
t0=1800;
L=（251-1）*dt;
t=1800:
2:
2300;
subplot（311）,
plot（t,y）;
xlabel（'时间/ms'）;ylabel（'振幅'）;
title（'原始数据图像'）;gridon;
y=fft（y,2*L）;%对输入信号的FFT
y0=fftshift（y）;%频谱中心化,matlab内部函数
M=length（y0）;%FFT频率轴点数（maybedifferent）
df=fs/（M-1）;%频率步长
f=（0:
M-1）*df;%频率采样序列（矢量）
f0=f-f（M/2+1）;
subplot（312）,plot（f0,y0）;%随频率变化的振幅
xlabel（'频率/Hz'）;ylabel（'振幅'）;
title（'对称频谱图：
'）;gridon;
xlim（[-fs/2,fs/2]）;
subplot（313）,plot（f（1:
M/2）,y（1:
M/2））;%绘制有效频谱
xlabel（'频率/Hz'）;ylabel（'振幅'）;
title（'只含正半轴频谱图：
'）;gridon;
xlim（[0,fs/2]）;
Butterworth）低通滤波器（一维）的频率响应曲线。
D0=100;
x=0:
50:
30000;
[M,N]=size（x'）;
N1=fix（M/2）;
n=[1,2,4];
y=zeros（3,M）;
h=zeros（1,M）;
forj=1:
3;
fori=1:
M
d=sqrt（（i-N1）^2）;
h（i）=1/（1+（d/D0）^（2*n（j）））;
end
y（j,:
）=h;
end
plot（x,y（1,:
））;
holdon
plot（x,y（2,:
））;
holdon
plot（x,y（3,:
））;
xlabel（'频率/HZ'）;ylabel（'增益'）;
title（'Butterworth低通滤波器（一维）的频率响应曲线'）;
gtext（'一阶曲线'）
gtext（'二阶曲线'）
gtext（'四阶曲线'）
2阶Butterworth低通滤波器，
clearall;
closeall;
I=importdata（'seismic.txt'）;
y=I（:
2）;
dt=2;
fs=1/0.002;
L=251;
t0=1800;
L=（251-1）*dt;
t=1800:
2:
2300;
subplot（411）,
plot（t,y）;
xlabel（'时间/ms'）;ylabel（'振幅'）;
title（'原始数据图像'）;gridon;
y=fft（y,2*L）;%对输入信号的FFT
y=fftshift（y）;
[M,N]=size（y）;
df=fs/（M-1）;%频率步长
f=（0:
M-1）*df;
subplot（412）
plot（f,y）;
xlabel（'频率/HZ'）;ylabel（'振幅'）;
title（'FFT变换所得频谱图'）;
D0=50;
n=2;
N1=fix（M/2）;
fori=1:
M
d=sqrt（（i-N1）^2）;
h=1/（1+（d/D0）^（2*n））;
y1（i,1）=h*y（i,1）;
end
subplot（413）
plot（f,y1）;
xlabel（'频率/HZ'）;ylabel（'振幅'）;
title（'Butterworth低通滤波所得曲线图1:
D0=50HZ'）;
D0=200;
fori=1:
M
d=sqrt（（i-N1）^2）;
h=1/（1+（d/D0）^（2*n））;
y1（i,1）=h*y（i,1）;
end
subplot（414）
plot（f,y1）;
xlabel（'频率/HZ'）;ylabel（'振幅'）;
title（'Butterworth低通滤波所得曲线图:
D0=200HZ'）;
四、思考题
1.分别阐述和解释什么叫信号的采样频率、奈奎斯特（Nyquist）频率、时间采样率及频率分辨率？
答：
采样频率，也称为采样速度或者采样率，定义了每秒从连续信号中提取并组成离散信号的采样个数。
奈奎斯特频率是离散信号系统采样频率的一半，因哈里·奈奎斯特或奈奎斯特－香农采样定理得名。
采样定理指出，只要离散系统的奈奎斯特频率高于被采样信号的最高频率或带宽，就可以真实的还原被测信号。
反之，会因为频谱混叠而不能真实还原被测信号
采样频率的倒数是采样周期或者叫作采样时间，它是采样之间的时间间隔。
通俗的讲采样频率是指计算机每秒钟采集多少个声音样本。
频率分辨率：
可以理解为在使用DFT时，在频率轴上的所能得到的最小频率间隔f0=fs/N=1/NTs=1/T,其中N为采样点数，fs为采样频率，Ts为采样间隔。
所以NTs就是采样前模拟信号的时间长度T，所以信号长度越长，频率分辨率越好。
2.根据所学知识，简要叙述离散傅立叶变换（DFT）和离散余弦变换（DCT）在数字信号处理中的主要用途。
答：
DFT是现代信号处理的基础，应用非常广泛，比如自适应滤波器啊，阵列信号处理、正交频分复用等等都用的到。
离散余弦变换（DCT）用于语音信号和图像信号的变化处理。
3.简要叙述频率滤波与时域滤波在处理上有什么不同。
答：
时域的卷积是在空间上图像矩阵进行卷积运算，这种过程，用一般的编程工具，就是双重嵌套循环，所以执行会比较耗时。
频域的话要先进行傅里叶变换，变换了之后，其实就是和相关的算子进行乘积运算，这个过程肯定比上述的卷积运算快很多。
但频域的处理较之时域的处理多了傅里叶和逆向傅里叶的计算，这部分开销同样也比较大的。
但频域的处理有时域不能比拟的优势，就是对于不规则的噪声值，通过傅里叶变换可以得到很好的平滑效果；但相应的，时域在边缘提取上，要比频域的处理更优秀。
4.实验内容1中，式
（1）表示的信号是平稳信号还是非平稳信号，简要说明平稳信号和非平稳信号的区别
答：
平稳信号
平稳信号即其一阶矩为常数，二阶矩与信号时间的起始点无关，只和起始时间差有关。
非平稳信号是指分布参数或者分布律随时间发生变化的信号
升级会员 