数字信号处理设计与仿真分析

本文的pdf文件：link

一、问题重述

1、建立两个模拟信号的数学模型Sa1(t)S_a1 (t)Sa1(t)和Sa2(t)S_a2 (t)Sa2(t)，其中Sa1(t)S_a1 (t)Sa1(t)为有用信号, Sa2(t)S_a2 (t)Sa2(t)为干扰信号。两个信号的中心频率，信号带宽等参数由学生自己选定，要求两个信号的频率不重叠，Sa2(t)S_a2 (t)Sa2(t)的幅度比Sa1(t)S_a1 (t)Sa1(t)的幅度高20dB，两个信号时域叠加得到合成信号Xa(t)X_a (t)Xa(t),即Xa(t)=Sa1(t)+Sa2(t)X_a (t)= S_a1 (t)+ S_a2 (t)Xa(t)=Sa1(t)+Sa2(t)，设计计算机程序仿真产生以上三个信号，分别画出三个模拟信号的时域波形和频谱图。
       2、根据Xa(t)X_a (t)Xa(t)的中心频率和带宽，按照奈奎斯特采样定理选择采样频率fs，分别对三个信号进行时域采样，得到离散信号S1(n),S2(n),x(n)S_1 (n), S_2 (n), x(n)S1(n),S2(n),x(n)。利用FFT算法分析离散信号的频谱，分别画出三个离散信号的时域波形和频谱图。
       3、设计数字滤波器H(z),要求该滤波器对干扰信号S2(n)S_2 (n)S2(n)的衰减大于40dB.提出滤波器的设计指标，并设计滤波器，给出滤波器的设计结果，绘制滤波器的幅频特性和相频特性曲线，验证滤波器的设计结果是否达到设计指标要求。
       4、选择实现数字滤波器H(z)的结构，画出结构信号流图
       5、将合成信号x(n)输入数字滤波器H(z),按照所选择的滤波器结构，设计计算机程序计算滤波器的输出响应y(n)，画出y(n)的时域波形和频谱图

二、相关术语

2.1 中心频率

周期性的信号均有其对应的频率，而且可以透过傅里叶级数转化为不同频率弦波的和。许多物理原件的特性会随着讯号的频率而改变，例如电容在低频时阻抗变大，高频时阻抗变小，而电感恰好相反，高频时阻抗变大，低频时阻抗变小。

2.2 信号带宽

信号带宽是信号频谱的宽度，也就是信号的最高频率分量与最低频率分量之差。信号带宽限定了允许通过该信道的信号的上限频率和下限频率，也就是限定了一个频率通带。

由信号频谱图可以观察到一个信号所包含的频率成分。把一个信号所包含谐波的最高频率与最低频率之差，即该信号所拥有的频率范围，定义为该信号的带宽。因此可以说，信号的频率变化范围越大，信号的带宽就越宽。

2.3 信号幅度

表示一个信号在不同频率下的幅值和相位，可以采用成为频谱图的表示方法。在傅里叶分析中，把各个分量的幅度|Fn|或 Cn 随着频率nω1的变化称为信号的幅度谱。

2.4 滤波器

滤波器是一种选频装置，可以使信号中特定的频率成分通过，而极大地衰减其他频率成分。利用滤波器的这种选频作用，可以滤除干扰噪声或进行频谱分析。换句话说，凡是可以使信号中特定的频率成分通过，而极大地衰减或抑制其他频率成分的装置或系统都称之为滤波器。

2.5 奈奎斯特频率

奈奎斯特频率（Nyquist frequency）是离散信号系统采样频率的一半，因哈里·奈奎斯特（Harry Nyquist）或奈奎斯特－香农采样定理得名。采样定理指出，只要离散系统的奈奎斯特频率高于采样信号的最高频率或带宽，就可以避免混叠现象。
从理论上说，即使奈奎斯特频率恰好大于信号带宽，也足以通过信号的采样重建原信号。但是，重建信号的过程需要以一个低通滤波器或者带通滤波器将在奈奎斯特频率之上的高频分量全部滤除，同时还要保证原信号中频率在奈奎斯特频率以下的分量不发生畸变，而这是不可能实现的。在实际应用中，为了保证抗混叠滤波器的性能，接近奈奎斯特频率的分量在采样和信号重建的过程中可能会发生畸变。因此信号带宽通常会略小于奈奎斯特频率，具体的情况要看所使用的滤波器的性能。

2.6 信号流图

借助拓扑图形求线性代数方程组解的一种方法。在1953年由S.J.梅森提出，故又称梅森图。这一方法能将各有关变量的因果关系在图中明显地表示出来，常用于分析线性系统。

三、问题分析

由问题一可知，我们首先拿到的是一个模拟信号，为了模拟实际中的信号，我们人为为这个信号加上了干扰，最后得到了合成信号，由于在时域中只能看到信号的波形大小，如果在时域对信号进行去噪处理，难度很大，于是将其转换到频域，在转换的过程中，为了方便计算机处理，我们首先将信号进行采样，然后对其进行傅里叶变换，发现噪声的频率高于真实信号的频率，于是我们想到了设计滤波器来实现高频噪声的滤除，根据噪声的各个参数，设计出合适的滤波器，最终得到了一个滤波器可以滤除类似的一系列噪声，并输出干净的有用信号。整个实验的设计流程如下所示：

四、问题求解

4.1 问题1

建立两个模拟信号的数学模型Sa1(t)S_a1 (t)Sa1(t) 和Sa2(t)S_a2 (t)Sa2(t)，其中Sa1(t)S_a1 (t)Sa1(t)为有用信号, Sa2(t)S_a2 (t)Sa2(t)为干扰信号。两个信号的中心频率，信号带宽等参数由学生自己选定，要求两个信号的频率不重叠，Sa2(t)S_a2 (t)Sa2(t)的幅度比Sa1(t)S_a1 (t)Sa1(t)的幅度高20dB，两个信号时域叠加得到合成信号Xa(t)X_a (t)Xa(t),即Xa(t)=Sa1(t)+Sa2(t)X_a (t)= S_a1 (t)+ S_a2 (t)Xa(t)=Sa1(t)+Sa2(t)，设计计算机程序仿真产生以上三个信号，分别画出三个模拟信号的时域波形和频谱图。

4.1.1 建立有用信号

首先生成一个有用的信号Sa1(t)S_a1 (t)Sa1(t)，并计算其傅里叶变换，查看其频谱，并根据其幅度，设计干扰信号，使其干扰信号的幅度高于有用信号的20dB,并使得两个信号的频率不重叠，最后将两个信号相加，得到合成信号。
       为了选择合适的函数，我们选择了双边指数函数作为原始信号，选择双边指数函数的原因为：其傅里叶变换后的频谱图大致呈现一个三角形，这样我们可以更好的选择干扰信号的频率。

       上式中α为系数，下面我们选择α=1，得到确定的函数，其图像如下图所示：

       对Sa1(t)S_a1 (t)Sa1(t)进行傅里叶变换，得到其频谱图，其傅里叶变换过程如下：

       带入α=1，得到其傅里叶变换后的函数，画出其图像如下：

4.1.2 建立干扰信号

因为干扰信号的幅度比有用信号的幅度高20dB，即

20log⁡(x)=20dB20 log⁡(x)=20dB 20log⁡(x)=20dB

求解得到x=10，有用信号的最大值为1，则干扰信号的最大值设为10，同时因为有用信号的频谱主要集中在0位置处，对干扰信号进行频谱的搬移，就可以使得两个信号的频谱不重叠。设置干扰信号为
Sa2(t)=cos⁡(ω0t)S_a2 (t)=cos⁡(ω_0 t)Sa2(t)=cos⁡(ω0t)

幅度扩大十倍得到

Sa2(t)=10cos⁡(ω0t)S_a2 (t)=10cos⁡(ω_0 t)Sa2(t)=10cos⁡(ω0t)

对上市进行傅里叶变换得到

F(jω)=10∗pi∗[δ(ω+ω0)+δ(ω−ω0)]F(jω)=10*pi*[δ(ω+ω_0 )+δ(ω-ω_0 )]F(jω)=10∗pi∗[δ(ω+ω0)+δ(ω−ω0)]

得到其频谱主要集中在ω0ω_0ω0处，为使得频谱进行搬移，对应的时域变化为

       取ω0=10ω_0=10ω0=10，得到最终的干扰信号为
       根据函数表达式，画出其时域波形如图所示：

       由图像可以看出，干扰信号的幅度为10，满足题目所给条件，对该信号进行傅里叶变换，得到其频域图如下：

       由图像可知，频域上该信号的主要频率集中在ω=10处，两个信号的频谱无混叠。

4.1.3 计算合成信号

根据已经得到的原始信号和干扰信号函数，可得合成信号为：

画出其时域图像如下：

对其进行傅里叶变换，得到其频谱图像如下：

4.2 问题2

根据Xa(t)X_a (t)Xa(t)的中心频率和带宽，按照奈奎斯特采样定理选择采样频率fs，分别对三个信号进行时域采样，得到离散信号S1(n),S2(n),x(n)S_1 (n), S_2 (n), x(n)S1(n),S2(n),x(n)。利用FFT算法分析离散信号的频谱，分别画出三个离散信号的时域波形和频谱图。

4.2.1 时域采样

根据信号的频率响应，计算奈奎斯特频率，采用该频率对模拟信号进行采样。

4.2.2 FFT变换

傅里叶变换的一般规律：

由于计算机无法有效的求解模拟信号的傅里叶变换，所以引入了离散时间傅里叶变换，首先将模拟信号进行时域采样得到离散时间信号，然后对采样信号进行DTFT，但是由于DTFT变换后得到的频域波形是连续的，依然不适合计算机处理，又引入了DFT，针对离散周期序列，就可以得到离散的频谱，随着FFT的出现，序列的傅里叶变化得到了极大的应用。

各类变换对应的频谱形式如下所示：

对得到的离散时间信号进行快速傅里叶变换，得到三个信号的频谱图。

4.3 问题3

设计数字滤波器H(z),要求该滤波器对干扰信号S2(n)S_2 (n)S2(n)的衰减大于40dB.提出滤波器的设计指标，并设计滤波器，给出滤波器的设计结果，绘制滤波器的幅频特性和相频特性曲线，验证滤波器的设计结果是否达到设计指标要求。

4.3.1 数字滤波器

指输入输出均为数字信号，通过一定运算关系改变输入信号所含频率成分的相对比例或者滤除某些频率成分的器件。根据特性可以分为低通，高通，带通和带阻滤波器。

       常用的数字滤波器一般属于选频滤波器，假设数字滤波器的频率响应为：

       其中，|H(e^jω )|称为幅频特性函数，表示信号通过该滤波器后各频率成分振幅衰减情况。Θ(ω)为相频特性函数，反映各频率成分通过滤波器后在时间上的延时情况。数字低通滤波器的幅频响应曲线如图：

       但是由于数字滤波器的设计不如模拟滤波器的技术成熟，所以IIR滤波器一般采用先将数字滤波器应有的参数转化为模拟滤波器响应的参数，之后根据模拟滤波器的模板设计模拟滤波器，在将其转换为数字滤波器，其设计流程如下：

4.3.2 matlab实现

Matlab信号处理工具箱中提供了设计巴特沃斯模拟滤波器的函数，各函数的基本格式如下：

根据matlab函数，认为设定滤波器的参数，得到合适的滤波器，得到其幅频响应和相频响应如下图所示：

根据设计的干扰信号和有用信号的特点，设计响应的滤波器，这里选择巴特沃斯滤波器，设计通带截止频率为3rad/s,阻带截止频率为7rad/s，通带最大衰减为1dB，阻带最小衰减为40dB，根据滤波器的幅频响应曲线，该滤波器满足设计要求。

4.4 问题4

选择实现数字滤波器H(z)的结构，画出结构信号流图。
       IIR数字滤波器的基本结构为：直接型结构，级联型结构，并联型结构。系统的差分方程形式如下：

       直接I型结构：

       直接II型结构：

       根据计算出来的参数，采用直接II型结构画出信号流图：
       根据求解得到滤波器的参数如下：
       N = 7
       Wc = 3.6257
       B = 1.0e+03 * 0 0 0 0 0 0 0 8.2358
       A =1.0e+04 * 0.0001 0.0016 0.0133 0.0695 0.2521 0.632 1.0208 0.8236
       根据求解得到的参数，可以画出信号流图，这里采用直接II型，如图所示：

4.5 问题5

将合成信号x(n)输入数字滤波器H(z),按照所选择的滤波器结构，设计计算机程序计算滤波器的输出响应y(n)，画出y(n)的时域波形和频谱图。
在设计好了模拟滤波器之后，我们需要将模拟滤波器转化为数字滤波器，主要有脉冲响应不变法和双线性变换法。

4.5.1 脉冲响应不变法

主要是时域逼近的思想，利用DF的单位脉冲响应h[k]模仿AF的单位冲激响应h(t)。

AF极点与DF极点的映射关系为：

根据此映射关系可以得到从S平面到Z平面的分布：

4.5.2 双线性变换法

基本思想是将非带限的模拟滤波器映射为最高角频率的带限模拟滤波器。

       利用脉冲响应不变法将上述模拟滤波器转换为数字滤波器。再将时域采样合成信号得到的离散信号输入上述滤波器，得到输出响应y(n)。
       其时域波形图，如下所示：

       由上述结果，可以发现：输入信号Xa(t)经过滤波器后，在生成的输出响应信号yn的频谱中，高频分量被明显滤除，只剩下了低频分量，说明所设计的滤波器符合低通要求。

五、结果分析

对比真实信号和经过滤波器输出的信号的频谱，我们可以发现，高频分量的频谱基本被滤除，最终得到了去噪之后的信号

本文中，选择的有用信号是一个双边指数序列，而他的傅里叶变换只一个类似三角形的分布，并且噪声的设置也是有意为之，很明显噪声的频率高于有用信号的频率，所有，简单的设置一个低通滤波器就可以将噪声滤除。但是现实生活中的信号形式更加多变，且噪声的频率并不一定明显高于有用信号，所以在设置滤波器时难度会更大，可能会经过很多次的滤波，才可以得到比较纯净的有用信号.

以上为常用的滤波器，在现实生活中的信号形式更加多变，且噪声的频率并不一定明显高于有用信号，所以在设置滤波器时难度会更大，应当根据实际情况选择合适的滤波器，还有可能会经过很多次的滤波，才可以得到比较纯净的有用信号。

六、matlab代码

syms t S1 S2;
% 有用信号
S1t=exp(-abs(t));
F1 = fourier(S1t);
figure(1);
%Sa1的时域波形和频谱图
fplot(S1t);
axis([-5,5,-0.5,1.2]);
title('S_a_1(t)的时域波形');
xlabel('t');
grid on;
figure(2)
fplot(F1);
title('S_a_1(t)频谱波形');
axis([-5,5,0,2.5]);
xlabel('ω');
grid on;% 干扰信号
S2t=10*exp(-abs(t)).*cos(10*t);
F2 = fourier(S2t);
%Sa2t的时域波形和频谱图
figure(3)
fplot(S2t);
axis([-3,3,-11,11]);
title('S_a_2(t)的时域波形');
xlabel('t');
grid on;
figure(4)
fplot(F2);
title('S_a_2(t)频域波形');
axis([-20,0,0,15]);
xlabel('ω');
grid on;%合成信号
Xat=S1t+S2t;
F3 = fourier(Xat);
%Xat的时域波形和频谱图
figure(5)
fplot(Xat);
axis([-2.5,2.5,-11,11]);
title('X_a(t)的时域波形');
xlabel('t');
grid on;
figure(6)
fplot(F3);
title('X_a(t)频域波形');
axis([-30,0,0,12]);
xlabel('ω');
grid on;%三个离散信号的时域波形和频谱图
N = 100;
n = 0:1:N-1;
fs = 30;
Ts = 1/fs;
S1n=exp(-abs(n*Ts));
S2n=10*exp(-abs(n*Ts)).*cos(10*n*Ts);
Xan=S1n+S2n;
S1k=abs(fft(S1n,N));
S2k=abs(fft(S2n,N));
Xak=abs(fft(Xan,N));
figure(7);
stem(n,S1n,'.');
grid on;
title('S_1(n)的时域波形');
xlabel('n');
ylabel('S_1(n)'); figure(8)
stem(n,S1k,'.');
title('S_1(n)的频谱图');
xlabel('k');
ylabel('S_1(k)');
grid on;figure(9)
stem(n,S2n,'.');
grid on;
title('S_2(n)的时域波形');
xlabel('n');
ylabel('S_2(n)');figure(10)
stem(n,S2k,'.');
title('S_2(n)频谱图');
xlabel('k');
ylabel('S_2(k)');
grid on;figure(11)
stem(n,Xan,'.');
grid on;
title('X_a(n)的时域波形');
xlabel('n');
ylabel('Xa(n)');figure(12)
stem(n,Xak,'.');
title('X_a(n)的频谱图');
xlabel('k');
ylabel('X_a(k)');
grid on;%由信号Xat的频谱所得
Wp = 3;
Ws = 7;
Rp = 1;
Rs = 40;
[N,Wc]=buttord(Wp,Ws,Rp,Rs,'s');
[B,A]=butter(N,Wc,'low','s');
[Bz,Az] = impinvar(B,A,fs);
W=0:0.001:50;
[H,W]=freqs(B,A,W);
phi=angle(H);
H=20*log10(abs(H)); figure(13);
plot(W,H);
title('滤波器的幅频响应');
xlabel('w');
ylabel('|Hw|');
grid on;figure(14)
plot(W,phi);
title('滤波器的相频响应');
xlabel('w');
ylabel('angle(Hw)');
grid on;yn = filter(Bz,Az,Xan);
n = 0:length(yn)-1;
Yk=abs(fft(yn,length(yn)));figure(16)
plot(n,yn);
title('滤波之后的y_n时域波形');
xlabel('n');
ylabel('yn');
grid on;figure(17)
stem(n,Yk,'.');
title('滤波之后的y_n的频谱');
xlabel('w');
ylabel('Yk');
grid on;