数字信号的 DFTDCT 及频域滤波.docx
《数字信号的 DFTDCT 及频域滤波.docx》由会员分享,可在线阅读,更多相关《数字信号的 DFTDCT 及频域滤波.docx(16页珍藏版)》请在冰豆网上搜索。
![数字信号的 DFTDCT 及频域滤波.docx](https://file1.bdocx.com/fileroot1/2022-12/16/9ee04878-753a-4b5f-8149-7b550a4ba3d4/9ee04878-753a-4b5f-8149-7b550a4ba3d41.gif)
数字信号的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)表示的信号是平稳信号还是非平稳信号,简要说明平稳信号和非平稳信号的区别
答:
平稳信号
平稳信号即其一阶矩为常数,二阶矩与信号时间的起始点无关,只和起始时间差有关。
非平稳信号是指分布参数或者分布律随时间发生变化的信号