数字信号处理实验(二)

实验目的

音频信号实验
产生两段音频信号，将声音信号进行傅立叶变换，将时域信号转为频域信号，并分别查看两段信号变化后的幅频特性，然后将两段频域信号线性叠加，合成新的频域信号，最后将此频域信号进行反傅立叶变换得到时域信号。
DFT实验
掌握DFT(FFT)对时域离散信号进行频谱分析的方法。
滤波器的设计
掌握窗口函数法设计FIR数字滤波器的原理和方法

实验原理

音频信号实验

利用matlab实现音频的读取，快速傅立叶变换，频域序列的叠加，傅立叶反变换。
核心函数为
wavrecord() 通过设备获取音频信号
wavread() 将音频文件读入
wavplay() 播放wav格式的音频时域信号
fft() 离散快速傅立叶变换
ifft() 离散快速傅立叶凡变换

2. DFT实验

DFT变换的公式为
$x(k)=\sum_{n=0}^N x(n)e^{-j kn{2\pi/N}}, k=0, 1, 2, ...N-1$ (1)
反变换的形式为
$x(n)=(1/N)\sum_{k=0}^N-1x(k)e^{jkn2\pi/N}, n=0, 1, 2, ..., N-1$(2)
Matlab 中计算DFT(FFT)的函数为fft() 反函数为ifft()
fft()函数的说明
通过查看说明文档，fft函数共有四种形式
Y = fft(x)
Y = fft(X, n)
Y = fft(X, [], dim)
Y = fft(X, n, dim)
Y = fft(x)返回向量x的离散傅立叶变换，如果x是个矩阵，那么返回每一行的离散傅立叶变换。
Y = fft(X, n) 返回n点的离散傅立叶变换。
ifft()函数的说明
y = ifft(X)
y = ifft(X, n)
y = ifft(X, [], dim)
y = ifft(X, n, dim)
y = ifft(…, ‘symmetric’)
y = ifft(…, ‘nosymmetric’)
y = ifft(X) 返回向量x的离散傅立叶反变换，如果x是个矩阵，那么返回每一行的离散傅立叶反变换。
y = ifft(X, n) 返回n点的离散傅立叶反变换。

3. 滤波器的设计

a. 窗函数法设计FIR数字滤波器
设欲设计的滤波器的理想频域相应为$H_d(e^{jw})$ ，单位脉冲响应为则$$h_d(n)=1/2\pi \int_{-\pi}^{\pi}H_d(e^{jw})dw$$根据给定的$H_d(e^{jw})$求出$h_d(n)$一般是无限长且非因果的。为了得到一个因果的有限长的滤波器$h(n)$ ，需要设计一个窗口函数$w(n)$对 $h_d(n)$进行截断处理
$h(n)=h_d(n)w(n)$

于是$h(n)$就是十几设计FIR滤波器的单位脉冲响应，其频域相应 $H(e^{jw})$为
$H(e^{jw})=\sum_{n=0}^{N-1}e^{-jwn}$
其中N为窗口$w(n)$的长度。窗口函数的形状和窗口长度N决定了窗口函数法设计出的FIR滤波器的性能。

实验内容及结果分析

音频信号实验

第一步，读入音频信号，并播放获取的信号代码如下

1.  [x1,FS1,NBITS1]=wavread('./sunhao1.wav'); %文件在当前文件夹
2.  wavplay(x1);
3.  pause(2);
4.  [x2,FS2,NBITS2]=wavread('./sunhao2.wav');
5.  wavplay(x2);
6.  pause(2);
7.  figure
8.  subplot(2,1,1),plot(x1);
9.  subplot(2,1,2),plot(x2);

两段音频文件的时域图谱为

第二步，声音信号进行离散傅立叶变换，对应代码如下

1.  y1=x1(:,1); y2=x2(:,1);
2.  N=min(length(y1),length(y2));
3.  z1=fft(y1,N); z2=fft(y2,N);
4.  k=0:N-1; wk=2*k/N;
5.  figure; subplot(2,1,1);
6.  stem(wk,abs(z1),'linewidth',2);
7.  title('(a) 音频1的幅频特性图');xlabel('omiga/pi');ylabel('幅度');grid
8.  subplot(2,1,2); stem(wk,abs(z2),'linewidth',2);
9.  title('(b) 音频2的幅频特性图');xlabel('omiga/pi');ylabel('幅度');grid

两段声音的幅频特性图谱

从图中可以看出在低频和高频处幅度较大，其他频率部分幅度较低，且相对平缓。
最后，将两段频域信号进行叠加，叠加后的序列再进行反傅立叶变换到时域信号。

DFT实验

第一部分
对矩形序列分别进行不同点的快速傅立叶变换，令 x=[1, 1, 1, 1]对x分别进行N=8，N=16, N=32的fft变换

结果如图所示，(a)和(b)子图表示的是N=8时的fft的幅度和相位特性图，(c)和(d)子图表示的是N=16时的fft的幅度和相位特性图；(e)和(f)子图表示的是N=32时的fft的幅度和相位特性图。可以看出N越大图形提供的信息越完整。
第二部分，对余弦序列进行快速傅立叶变换。

图中(a)和(b)子图表示的是N=8时的fft的幅度和相位特性图，(c)和(d)子图表示的是N=16时的fft的幅度和相位特性图；(e)和(f)子图表示的是N=32时的fft的幅度和相位特性图。虽然看不出明显的差别。

滤波器的设计

窗函数法设计线性相位FIR滤波器
(1) 数字滤波器的性能要求：临界频率$w_0$，滤波器单位脉冲响应长度N
(2) 根据性能要求，合理选择单位脉冲响应$h(n)$的奇偶对称性，从而确定理想频率响应 $H_d(e^{jw})$的幅频特性和相频特性。
(3)求理想单位脉冲响应$h_d(n)$ 。在实际中，可对 $H_d(e^{jw})$按M(M远远大于N)点等距采样，并对其求IDFT得 $h_M(n)$，用 $h_M(n)$代替 $h_d(n)$。
(4) 选择适当的窗函数$w(n)$ ，根据 $h(n)=h_d(n)w(n)$求所需设计的FIR滤波器单位脉冲响应。
(5) 求 $H(e^{jw})$，分析其幅频特性，若不满足要求，可适当改变窗函数形式或长度N，重复上述操作，直至得到满意结果。
常见的窗函数有以下几种形式

窗函数	旁瓣峰值幅度/dB	过度带宽	阻带最小衰减/dB
矩形窗	-13	4π/ N	-21
三角形窗	-25	8π/ N	-25
汉宁窗	-31	8π/ N	-44
哈明窗	-41	8π/ N	-53
布莱克曼窗	-57	12π/ N	-74
凯塞窗(a=7.865)	-57	10π/ N	-80

用窗口函数设计数字滤波器。
函数fir1()

1.  b = fir1(n, Wn)
2.  b = fir1(n, Wn, ‘ftype’)
3.  b = fir1(n, Wn, window)
4.  b = fir1(n, Wn, ‘ftype’, window)
5.  b = fir1(… , ‘nomalization’)

b = fir1(n, Wn)返回6dB截止频率为Wn的n阶（单位脉冲响应h(n)长度N=M+1）FIR低通（Wn为标量）滤波器系数向量b，默认选用哈明窗。当Wn是一个向量时，滤波器是带通滤波器。
fir1() 应用传统的线性相位方法设计滤波器。它可以设计标准的低通滤波器，高通滤波器，带通滤波器和傣族滤波器。该函数默认的初始化的滤波器的回馈大小中心值是0dB。
ftype有如下几种形式
‘high’ 高通滤波器具有Wn的截频.
‘stop’ 带阻滤波器，如果Wn=[w1, w2], 带阻的频率值被限制在区间之间。
‘DC-1’ 为多带滤波器设置第一个带。
‘DC-0’ 为多带滤波器设计阻带。
设计低通滤波器的演示程序如下：

实验拓展

设计低通滤波器，用正弦曲线 $f(t)=2sin(2\pi 20t) + 4sin(2\pi 60t)$来测试滤波效果。
代码如下

1.  t = 0:0.01:2;
2.  f =2*sin(2*pi*20*t)+4*sin(2*pi*60*t);  %原始图像
3.  N = 11;                %滤波器节点个数
4.  wc = 0.5;              %归一化截止频率
5.  hd = fir1(N,wc,'low'); % 基于加窗函数的FIR滤波器设计
6.  ft = conv(f,hd);
7.  figure(1)
8.  plot(abs(fft(f)));
9.  title('原始信号f');
10. figure(2)
11. plot(abs(fft(ft)));
12. title('滤波后信号ft');

如图所示(a)图表示原始图像，(b)图表示低通滤波后的图像，原始图像的函数原型为
。可以看到50-150之间两段频率较高的峰值部分被滤波器滤掉，剩下低频部分。

总结

在学习方法上，我有这点体会：学习工科，重在物理意义的理解。对于任何知识点，首先要尝试去理解这个知识点所表达的物理意义是什么，不要一开始就掉进了数学推导的茫茫大海中。先抓住主干，再去考量细节分支，最后再补充特殊情况。这是学习一个已经较为系统的知识的比较好的方法。若一开始从各种细节做起，则会茫然无头绪。
针对数字信号处理这门课程（目前只看到了DFT, FFT,后面的各种滤波器神马的还没有看。所以只拿DFT,FFT 说事儿。），我认为主干是这样的：每个信号都有一个频域特性，我们可以使用各种数学方法来观察信号的频域特性，不同的数学方法观察到的频域特性可能有所不同。这些数学方法包括：傅里叶变换（FT），离散时间傅里叶变换(IDFT)，离散傅里叶变换(DFT)和快速傅里叶变换(FFT)。在这四中数学处理方法中，只有DFT 和FFT 是可以在计算机中处理的，因为DFT和FFT是数值化的计算方法，而FT 和 IDFT是积分的计算方法。

附录

下面是实验过程中的代码

%#### 音频信号实验
%Step1 输入音频信号
%####
[x1,FS1,NBITS1]=wavread('./sunhao1.wav');
wavplay(x1);
pause(2);
[x2,FS2,NBITS2]=wavread('./sunhao2.wav');
wavplay(x2);
pause(2);
figure
subplot(2,1,1),plot(x1);
subplot(2,1,2),plot(x2);
%####
%Step2 输入音频信号
%####
y1=x1(:,1);
y2=x2(:,1);
N=min(length(y1),length(y2));
z1=fft(y1,N);
z2=fft(y2,N);k=0:N-1;
wk=2*k/N;
figure;
subplot(2,1,1);
stem(wk,abs(z1),'linewidth',2);
title('(a) 音频1的幅频特性图');xlabel('omiga/pi');ylabel('幅度');grid
subplot(2,1,2);
stem(wk,abs(z2),'linewidth',2);
title('(b) 音频2的幅频特性图');xlabel('omiga/pi');ylabel('幅度');gridz3=z1+0.3*z2;
z4=ifft(z3);
wavplay(z4);

%FIR 实验
xn = [1 1 1 1];
X_8 = fft(xn, 8);
X_16 = fft(xn, 16);
X_32 = fft(xn, 32);
figure('name', '矩形序列的fft频谱');k = 0 : 7;
wk = k * 2 / 8;
subplot(3, 2, 1); stem(wk, abs(X_8), 'LineWidth', 2);
title('（a）8点DFT的幅频特性图'); xlabel('oumege / pi'); ylabel('幅度'); grid on;
subplot(3, 2, 2); stem(wk, angle(X_8), 'LineWidth', 2);
title('(b) 8点DFT的相频特性图');grid on; xlabel('oumege / pi');ylabel('相位'); axis([0, 2, -3.5, 3.5]);k = 0 : 15; wk = 2 * k / 16;
subplot(3, 2, 3); stem(wk, abs(X_16), 'LineWidth', 2);
title('(c) 16点的幅频特性图'); xlabel('oumege / pi'); ylabel('幅度'); grid on;
subplot(3,2, 4); stem(wk, angle(X_16), 'LineWidth', 2);
title('(d) 16点的相频特性图'); grid on; xlabel('oumege / pi'); ylabel('相位');
axis([0, 2, -3.5, 3.5]);k = 0 : 31; wk = 2 * k / 32;
subplot(3, 2, 5); stem(wk, abs(X_32), 'LineWidth', 2);
title('(c) 32点的幅频特性图'); xlabel('oumege / pi'); ylabel('幅度'); grid on;
subplot(3,2, 6); stem(wk, angle(X_32), 'LineWidth', 2);
title('(d) 32点的相频特性图'); grid on; xlabel('oumege / pi'); ylabel('相位');
axis([0, 2, -3.5, 3.5]);clear all; n = 0 : 50;
xn = cos(pi / 4 * n);
X_4 = fft(xn, 4);
X_8 = fft(xn, 8);
X_16 = fft(xn, 16);
figure('name', '余弦序列的频谱');k = 0 : 3; wk = 2 * k / 4;
subplot(3, 2, 1); stem(wk, abs(X_4), 'LineWidth', 2);
title('(c) 4点的幅频特性图'); xlabel('oumege / pi'); ylabel('幅度'); grid on;
subplot(3,2, 2); stem(wk, angle(X_4), 'LineWidth', 2);
title('(d) 4点的相频特性图'); grid on; xlabel('oumege / pi'); ylabel('相位');
axis([0, 2, -3.5, 3.5]);k = 0 : 7; wk = 2 * k / 8;
subplot(3, 2, 3); stem(wk, abs(X_8), 'LineWidth', 2);
title('(c) 8点的幅频特性图'); xlabel('oumege / pi'); ylabel('幅度'); grid on;
subplot(3,2, 4); stem(wk, angle(X_8), 'LineWidth', 2);
title('(d) 8点的相频特性图'); grid on; xlabel('oumege / pi'); ylabel('相位');
axis([0, 2, -3.5, 3.5]);k = 0 : 15; wk = 2 * k / 16;
subplot(3, 2, 5); stem(wk, abs(X_16), 'LineWidth', 2);
title('(c) 16点的幅频特性图'); xlabel('oumege / pi'); ylabel('幅度'); grid on;
subplot(3,2, 6); stem(wk, angle(X_16), 'LineWidth', 2);
title('(d) 16点的相频特性图'); grid on; xlabel('oumege / pi'); ylabel('相位');
axis([0, 2, -3.5, 3.5]);

%拓展实验
t = 0:0.01:2;
f =2*sin(2*pi*20*t)+4*sin(2*pi*60*t);
N = 11;                %滤波器节点个数
wc = 0.5;              %归一化截止频率
hd = fir1(N,wc,'low'); % 基于加窗函数的FIR滤波器设计
ft = conv(f,hd);
figure(1)
subplot(1, 2, 1); plot(abs(fft(f)));
%title('原始信号f');
%figure(2)
subplot(1, 2, 2); plot(abs(fft(ft)));
%title('滤波后信号ft');