实验三DFT和DCT及频域滤波.docx

1、实验三DFT和DCT及频域滤波一 实验名称：数字信号的 DFT/DCT 及频域滤波二 实验目的1 熟练掌握数字信号（1D）及数字图像(2D)离散傅立叶变换（DFT）及离散余弦变换（DCT）方法、基本原理及实现流程。熟悉两种变换的性质，并能对 DFT 及 DCT 的结果进行必要解释。2 深入理解离散信号采样频率、奈奎斯特频率及频率分辨率等基本概念，弄清它们之间的相互关系。了解离散傅里叶变换（DFT）中频率泄露的原因，以及如何尽量减少频率泄露影响的途径。3 熟悉和掌握利用 MATLAB 工具进行 1D/2D FFT 及 DCT 的基本步骤、MATLAB 函数使用及对具体变换的处理流程。4 能熟练应

2、用 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，对信号进行离散化，并

3、画出一个周内的信号振幅随时间变化的波形图。2) 对离散信号进行傅立叶变换，分别画出频谱中心化及有效频率范围（不含负频）2种方式下的幅值()随频率 变化的分布图，要求纵横坐标正确标注物理量和单位。3) 对式（1）信号，加随机噪声，重复步骤（1）和（2）的处理过程。4) 通过对变换结果的分析，说明采样频率、奈奎斯特（Nyquist）频率（）及采样时间间隔T 三者之间的相互关系，并简要描述模拟信号的采样定理。图 1注：加随机噪声的处理只需改变初始信号即可，其它步骤同理。2. 数字音频信号的DFT1) 读取一段 0.5s 的预先录制的数字音频信号（ “yes.wav”或“no.wav”文件中任选其一）

4、，画出随时间变化的声波波形图。2) 对数字音频信号进行离散傅立叶变换（DFT） ，分别画出频谱中心化及有效频率范围(不含负频)2 种方式下的幅值（）随频率变化的分布图，要求纵横坐标正确标注物理量和单位。图 23. 数字音频信号的 DCT 和 IDCT1) 对上述音频信号做离散余弦变换（DCT），画出 DCT 变换系数（变换结果）图，并对变换结果进行必要的解释，说明 DCT 变换的主要用途。要求按 DCT 原理自行编写实现代码，不允许直接调用 MATLAB 的 dct()函数。2) 按原理自行一段 MATLAB 代码，对第（1）步处理结果进行离散余弦反变换(IDCT)，将计算结果与原始音频信号进

5、行比较，检验编写代码的正确性。3) 编写一段 MATLAB 代码，利用快速傅立叶变换（FFT）程序实现快速 DCT 算法（FCT），并将计算结果与直接调用 dct（）的处理结果进行比较，检验编写代码的正确性。 图 34. 综合应用题：实际信号的频谱分析及频域滤波1) 编写一从保存在本地磁盘的文本文件中读入一实际数字信号，已知该信号的时间采样率为 dt = 2ms。文件中的信号由 301个等长的按列排列的一维列信号组成，每个一维列信号有 251 个采样点，信号实际计时起点为 1800ms，延时长度为 L= (251-1)* 2ms =500ms。请读出其中的某一列信号，并画出该信号振幅随时间变化

6、的波形图，以 ms 为时间单位。2) 对第一步中抽取的其中一列信号做快速傅里叶变换（FFT），分别画出频谱中心化的对称频谱和只含有正半抽的信号频谱图，并对该信号做简要的频谱分析。要求规范的标注纵横坐标实际物理量和对应的单位。3) 设定截止频率 D0=100，试在同一张图上以不同线型画出 n =1，2，4 阶下的巴特沃思（ Butterworth）低通滤波器（一维）的频率响应曲线。要求标注规范地纵横坐标实际物理量和对应的单位。4) 选择合适的 D0，利用上述 2 阶 Butterworth 低通滤波器，对第（1）步读取的列信号进行滤波实验。并分析截止频率对滤波效果的影响。图 4五 实验结果及分析

7、1. 图 51D数字信号的FFT及频谱分析 图 6分析：通过对变换结果的分析采样频率、奈奎斯特频率及采样时间间隔三者之间有如下的相互关系 模拟信号的采样定理：当采样频率大于连续信号最高频率时，采样之后的数字信号完整地保留了原始信号中的信息，是连续信号离散化的基本依据。2. 图 7数字音频信号的DFT 分析：从上面的频谱可以看出该音频的频率成分主要集中在700Hz以内，共有4支峰。3. 图 8数字音频信号的 DCT 和 IDCT 分析：DCT变换将信号从时域转换到变换域上，通过对变换后的系数分析，原能量集中在少数系数上，可以提高编码效率，压缩数据。4. 图 9综合应用题：实际信号的频谱分析及频域

8、滤波图 10图 11图 12 该信号的第二列数据的有效频率主要集中在0-100Hz以内 通过设置不同截止频率的Butterworth低通滤波器，可以看出对于低频信号，截止频率从100Hz减小，原信号的有效频率则被滤去的变多，若是图像信号，则表现出图像变的模糊。六.实验心得体会和建议 心得体会：通过这次实验使我深刻了解了奈奎斯特定理、DFT和DCT的基本原理以及巴特沃斯低通滤波器的构造，与此同时在上机实验中熟悉了MATLAB编写FFT和DCT 的基本方法及步骤。在对相关频谱的分析过程中更加深刻的理解了原理，及相关用途。 建议：可以让大家不调用FFT函数，直接编写DCT。七.程序源代码1. 1D数

9、字信号的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频率轴点数(may be different)f = (0:M-1)*fs/M; % 频率采样序列（矢量）%=fchar = num2str(fs); % 采样率转

10、化为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);grid on;f0 = f-f(M/2);subplot(312),plot(f0,(mag0)*2/N); % 随频率变化的振幅xla

11、bel(频率/Hz);ylabel(振幅);title(全部频率：,ltext);grid on;subplot(313),plot(f(1:M/2),(mag(1:M/2)*2/N); % 绘制有效频谱xlabel(频率/Hz);ylabel(振幅);title(有效频率：,ltext);grid on;注：对于加性噪声只需将x变为x=2*sin(30*pi*t)+0.5*cos(120*pi*t)+4*sin(240*pi*t)+randn(size(t);即可，其它程序同理。 2. 数字音频信号的DFTx,fs=wavread(yes.wav);N=4000;n=0:N-1;t=n/fs

12、;plot(t,x),%画出一个周期内的信号振幅随时间变化的波形图%=y=fft(x);mag=abs(y);y=fftshift(y);%频谱中心化mag0=abs(y);%=M = length(y); % fft频率轴点数(may be different)f = (0:M-1)*fs/M; % 频率采样序列（矢量）%=fchar = num2str(fs); % 采样率转化为char(本文)nchar = num2str(N); % 样点数转化为char(本文)ltext = strcat(fs=,fchar,Hz,N=, nchar, points);%拼title字符%=subpl

13、ot(311),plot(t,x); % 随时间变化的振幅xlabel(t/s);ylabel(振幅);xlim(0,0.5);title(yes);grid on;f0 = f-f(M/2);subplot(312),plot(f0,(mag0)*2/N); % 随频率变化的振幅xlabel(频率/Hz);ylabel(振幅);title(全部频率：,ltext);grid on;subplot(313),plot(f(1:M/2),(mag(1:M/2)*2/N); % 绘制有效频谱xlabel(频率/Hz);ylabel(振幅);title(有效频率：,ltext);grid on;3.

14、 数字音频信号的 DCT 和 IDCTy,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);for x=1:N; for u=2:N; G(x,u)=cos(2*(x-1)+1)*(u-1)*pi/(2*N); endendG(1,:)=sqrt(2/N); G(1,1)=sqrt(1/N);F=G*y;subplot(312),plot(1:N)/Fs,F)title(DCT coeffcients);xla

15、bel(t/s);ylabel(振幅);f=G*F;subplot(313),plot(t,f/Fs*4);title(The recovery of the signal from DCT coefficients)xlabel(t/s);ylabel(振幅);4. 综合应用题：实际信号的频谱分析及频域滤波clc,clear,close all;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=t

16、0+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),列);grid on;%=y=fft(sl,p-1);mag=abs(y);y=fftshift(y);%频谱中心化mag0=abs(y);%=M= length(y); % fft频率轴点数(may be different)f = (0:M-1)*fs/

17、M; % 频率采样序列（矢量）%=%画出中心频率图f0 = f-f(M/2);figure;subplot(211),plot(f0,mag0); % 随频率变化的振幅xlabel(频率/Hz);ylabel(振幅);title(全部频率);grid on;%=subplot(212),plot(f(1:M/2),mag(1:M/2); % 绘制有效频谱xlabel(频率/Hz);ylabel(振幅);title(有效频率);grid on;%=%生成Butterworth矩阵H=zeros(3,M);JS=1 2 4;D0=100;for n=1:M H(1,n)=1/(1+(n-(M)/2

18、)/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);hold on;subplot(312),plot(f0,abs(H(2,:),-b);hold on;subplot(313),plot(f0,abs(H(3,:),-b);hold

19、on;xlim(f0(1),f0(M);xlabel(f0/Hz);%=z1=y.*H(2,:);figure(4),plot(f0,abs(z1),grid on;八.思考题1. 分别阐述和解释什么叫信号的采样频率、奈奎斯特（Nyquist）频率、时间采样率及频率分辨率？a) 信号的采样频率：对模拟信号进行A/D采样时，每秒钟对信号采样的点数。b) 奈奎斯特频率：指最低允许的抽样率，是带限信号频率宽度的2倍。c) 采样时间： 每进行一次采样所需要的时间。d) 频率分辨率：指将两个相邻谱峰分开的能力，也指分辨两个不同频率信号的最小间隔。2. 根据所学知识，简要叙述离散傅立叶变换（DFT）和离散

20、余弦变换（DCT）在数字信号处理中的主要用途。主要用途：a) DFT：频域滤波，图像增强，频谱分析。b) DCT：数据压缩，提高编码效率。3. 简要叙述频率滤波与时域滤波在处理上有什么不同。区别：对于不规则的噪声值，通过FFT进行频率滤波可以得到很好的平滑效果，但是FFT不可避免的会存在频谱泄露；而时域在边缘提取上，要比频域的处理更优秀。4. 实验内容 1 中，式（1）表示的信号是平稳信号还是非平稳信号，简要说明平稳信号和非平稳信号的区别。a) 式（1）表示的信号是平稳信号b) 区别：通过对统计特征的对比，非平稳信号的统计特性（均值、方差等）随着时间变化而变化，而平稳信号的统计特性不随时间变化。报告评分： 指导教师签字：