信号时频域分析和滤波

项目名称：信号时频域分析和滤波

源文件和代码数字信号处理项目报告-CSDN博客

1. 项目目标

1.1.1知识目标

1）掌握信号的时域和频域采样定理；

2）掌握信号的频域分析理论；

3）掌握数字滤波器的设计方法；

4）掌握信号的滤波方法。

1.1.2能力目标

1）能够运用数字信号的时频域分析和快速计算的方法，根据要求设计合理的数字信号处理相关算法，并合理选择参数；

2）能够根据实验数据的情况设计实验方案，正确地采集测试数据，构建测试处理方法开展测试。

3）能够构建数字滤波器，对获取的实验数据进行时域和频域分析解释，进一步研究滤波器参数的选择和优化问题，提高信号处理算法的可靠性和有效性。

1.1.3素质目标

1）能够与小组同学合作；能够撰写报告和在答辩中清晰阐明自己的观点；

2）能够在实践过程中体会数字信号处理方法的精妙之处，每一个创新方法的起源，并不是离我们很远，创新方法的提出对社会的进步、国家的强盛起到了怎样的作用，激发学生的学习兴趣、对创新的渴望、使命担当的爱国精神。

2. 项目内容与要求

本项目要求完成信号采集，并对所采集信号进行时域、频域分析以及滤波处理。在报告中要求给出所实现功能的理论依据及具体实现方法，给出相应程序（附程）、中间处理结果和最终结果， 并对得到的结果进行分析，得出结论。

2.1 一维离散时间信号采集和分析

2.1.1 信号的时域分析

（1）采集离散一维时间信号

一维连续时间信号首先需要对其进行采样：

1）矩形脉冲信号、余弦信号；

2）实际的信号，例如语音信号等。

（2）画出所采集到的一维信号的时域波形

绘制出采集信号的时域波形，采用不同的采样间隔，注意采样间隔是时间。

图1  16HZ采样20点的余弦信号

图2  16HZ采样15点的余弦信号

离散抽样间隔大小的选择需要根据被采样信号的频率进行决定。根据奈奎斯特采样定理，为了避免采样信号时出现混叠现象，采样频率应该大于等于被采样信号的最高频率的两倍。因此，抽样间隔大小应该小于等于被采样信号周期的一半。根据采样定理，如果信号是带限的并且采样频率高于信号带宽的两倍，那么原来的连续信号可以从采样样本中完全重建出来。

2.1.2 信号的频域分析

（1）计算并画出信号的DTFT幅值谱结果波形

1）理论依据

DTFT的计算公式为

                         （1）

2）实现提示

1. 数字角频率

   

   的变化是连续的， 数字角频率的变化是离散的，离散间隔的选择取决于采样频率和采样点数。在DFT运算中，离散间隔

为2π/N，其中N为采样点数。

离散时间信号的频域抽样定理，对于一个离散时间信号，如果它的采样频率大于等于信号最高频率的两倍，那么就可以通过这个采样信号的频域信息完全还原原始信号。栅栏效应是指，在进行离散时间信号的频域采样时，由于采样点的数量有限，会在频域上形成一些虚假的频率分量，这些分量被称为栅栏效应。这会导致频域信息的失真和信号重构的误差。解决方法：增加采样点的数量，这样可以提高频率分辨率，减小栅栏效应；使用窗函数对信号进行加窗处理，这样可以减小频域泄漏，降低栅栏效应；对信号进行插值处理，这样可以增加采样点的数量，减小栅栏效应。

（2）要实现公式（1）的计算，观察需要的变量

x为一维时间信号

，

，w表示数字角频率

，dw表示数字角频率的步长。

DTFT的特点是它是一种连续的傅里叶变换，适用于离散时间信号的频谱分析。频谱范围从-π到π。

注意：离散时间信号频谱是以2π为周期的，我们只要计算一个周期即可， 为了便于观察，范围可选择两种情况，0到2π 或-π到π。

（3）程序实现提示：

N=length(x); （输入的序列的长度）

n=0:N-1;

dw=0.01*pi;（表示角频率的间隔，可选不用的间隔进行观察，也可以根据频域抽样定理来选dw=2*pi/N）

w=-pi:dw:pi （也可以取0到2π，观察与取-π到π的区别，为什么？这里要和后面DFT统一起来）

x，n，w都是一个1*N的行向量，这样（1）可由下面的这一个语句完成。

Xf=x*exp(-j*n'*w);

   答：n'为n的转置，转置后为列向量，与w行向量相乘后为N*N的矩阵，而x也为行向量，与exp(-j*n'*w)N*N矩阵相乘后为1*N行向量。

（2）计算并画出不同点数的DFT结果波形

1）理论依据

有限长序列

 

的离散时间傅里叶变换

 

在频率区间

 

的N个等间隔分布的点

 

上的取样值可以由下式表示：

  （2）

由上式可知，序列

 

的N点DFT

 

实际上就是

 

序列的DTFT在N个等间隔频率点

 

的抽样值

，所以可以利用DFT计算近似得到DTFT。

但当序列

 

的长度N较小时，在单位圆上的采样点比较稀疏（间距为2π/N），不能很好地表征其频响特性（栅栏效应），改进办法？

改进办法就是增加抽样点数N，也就是增加序列的长度，主要途径：

① 增加采样点，增加数据长度；

这种情况改变了原始序列，也就是对应长序列的情况，就不存在近似不足的情况，所以重点采用下面方法解决问题。

② 补零增加数据长度

这种情况未改变原始序列

，只是补零增加序列的长度，使得在单位圆上的采样点增多，得到的DFT谱线就更加精细，其包络就越接近DTFT的结果，这样就可以利用DFT计算DTFT。

2）实现提示

由于FFT是DFT的快速算法，因此在程序中可以直接使用MATLAB内置函数fft进行一维信号的DFT计算。

x=[2 -1 1 1];

Xk=fft(x); %计算序列x的DFT，序列长度是多少就计算了x几点的DFT

x=[2 -1 1 1 zeros(1,60)];%通过补60个零计算了x的64点DFT

Xk=fft(x);%注意这里计算的结果是复数

%也可直接使用：

Xk=fft(x,64);%(合并前面两句)

注意：

这里计算的DFT是在单位圆上的0-N-1的N个采样点，即对应0，

,…,

处的值。

计算结果是复数幅值谱和相位谱可以用fft快速傅里叶变换来求。

3）具体要求

计算不同点数DFT的结果，并画出幅值波形，与（1）得到的结果进行比较并分析得出结论。

当采样点数越多时，DFT的采样结果和DTFT的结果月接近。

图3     21点DFT与DTFT的幅值谱

  

图4      41点DFT与DTFT的幅值谱

2.2 加噪信号的时频域分析

2.2.1 噪声信号的时频域分析

主要选用下面两类噪声（加性）进行分析：

（1）正弦信号噪声（谐波信号）

（2）高斯白噪声

注意：不限于这两类噪声。

具体要求：

对这两类噪声信号进行介绍，画出噪声信号的时域和频域波形，观察分析噪声信号的时域和频域特点。

图5      正弦信号噪声的时域波形

图6      正弦信号的频域波形

图7       高斯白噪声时域波形

图8       高斯白噪声频域波形

高斯白噪声的时域特点为均值为0，方差为常数，频域中功率谱密度在整个频域内均匀分布。正弦信号的时域特点为振幅、频率、初相位和周期，频域中只有一个频率分量。

2.2.2 加噪信号的时域和频域分析

分别选一段规则信号（例如正弦信号）和语音信号添加噪声（加性），选择的信号要尽可能有代表性，长度合理；画出添加噪声后信号的时域、频域波形，分析发生了什么样的变化，语音主观听起来怎样，客观情况下信噪比。

当snr=3.6544时，可以听到类似耳鸣的声音，但语音信号的声音感觉更高，没有模糊的感觉。当snr=-4.8650时，可以明显听到类似耳鸣的声音，感觉已经高出了语音信号，但是语音信号依旧可以听清。

信噪比定义为接收到的无线信号与噪底之间的差值

实现提示：

y=x+noise1;%把噪声信号加入原始信号，noise1为噪声信号

正弦信号信噪比SNR=10，频率f=5HZ

图9    正弦信号加噪声后的频域波形

图10  正弦信号与正弦信号加噪声后的时域波形

图11    语音信号的时域波形

图12     语音信号加噪声后的时域波形（SNR=3.6544）

图13   语音信号加噪声后的时域波形（SNR= -4.8650）

图14     语音信号加噪声后的频域波形（SNR=3.6544）

图15     语音信号加噪声后的频域波形（SNR= -4.8650）

2.3 数字滤波器设计与实现

2.3.1 IIR滤波器设计方法和实现

IIR滤波器的设计选用巴特沃斯滤波器进行设计，可以选择冲激响应不变法或双线性变换法进行模拟滤波器到数字滤波器的转换。

（1）理论依据

1）冲激响应不变法

就是使数字滤波器的冲激响应序列

等于模拟滤波器的冲激响应

的采样值，即

             （3）

式中，T为采样周期。

已知

求

的方法是：

•  求模拟滤波器的冲激响应

  

  ：

                      （4）

•  求数字滤波器的冲激响应序列

  

                  （5）

•  求数字滤波器的系统函数：

                 （6）

2）双线性变换法

由于s平面和z平面的单值双线性映射关系为：

；   

                  （7）

式中，T为采样周期。

则已知

求

的方法是：

                 （8）

（2）IIR滤波器设计实现提示（巴特沃斯）

Wp=1/pi;Ws=2/pi; %通带截止频率、阻带截止频率；（归一化）根据实际情况给出。

f_p=5000; f_s=6500; Wp=f_p/(fs/2);Ws=f_s/(fs/2); %也可以根据模拟域情况进行换算，但须注意这是对应冲激响应不变法的情况，双线性变换需要采用不同的公式；

Rp=3; Rs=25; %通带波纹指数（通带最大衰减）、阻带最小衰减；（这要看加噪信号的实际情况进行设置）

[N,wc]=buttord(Wp,Ws,Rp,Rs);%巴特沃斯滤波器设计，计算滤波器阶数N和3dB截止频率wc。

[b,a]=butter(N,wc,'low');%计算滤波器系统函数分子分母多项式系数

freqz(b,a,1000,fs);%绘制滤波器特性曲线

（3）滤波器设计实验

1）取

Wp=0.1*pi;Ws=0.2*pi; Rp=3; Rs=25; fs=15000;

经计算所设计的低通滤波器：

N=4.2，取N=5；Wc=0.3142。

使用MATLAB计算得到的结果如图：

N=5；Wc=0.3534。与计算结果无明显差异

2）由1可知设计的滤波器满足性能要求。

使用双线性法得到的结果如图：即令

Wp=（2/T）*tan(Wp/2); Ws=（2/T）*tan(Ws/2);

3）它们的主要区别在于数字滤波器的频率响应与模拟滤波器的频率响应之间的映射方式不同。冲激响应不变法是将模拟滤波器的冲激响应直接映射到数字滤波器的冲激响应上。由于数字滤波器的采样频率有限，会导致模拟滤波器的高频部分被截断，从而影响数字滤波器的频率响应。

双线性变换法是将模拟滤波器的传递函数映射到数字滤波器的传递函数上。由于双线性变换是非线性的，会导致数字滤波器的幅值谱和相位谱与模拟滤波器的不同。

具体来说，双线性变换会将模拟滤波器的频率响应进行拉伸和压缩，从而导致数字滤波器的幅值谱和相位谱与模拟滤波器的不同。这种差异主要是由于数字滤波器的采样频率有限和双线性变换的非线性特性所导致的。

图16     模拟低通滤波器幅频特性

图17     双线性法低通滤波器幅频特性曲线

2.3.2 窗函数法设计FIR滤波器

利用窗函数法设计能够滤出噪声的FIR低通滤波器。

（1）理论依据

如果期望滤波器理想的频率响应函数为

，则其对应的单位脉冲响应为

             （9）

窗函数设计法的基本原理是用有限长单位脉冲响应序列

逼近

。由于

往往是无限长序列，而且是非因果的，所以用窗函数

将

截断，并进行加权处理，得到：

                    （10）

就作为实际设计的FIR数字滤波器的单位脉冲响应序列，其频率响应函数

为

                （11）

式(12)中，N为所选窗函数

的长度。

用窗函数法设计的滤波器性能取决于窗函数

的类型及窗口长度N的取值。设计过程中，要根据对阻带最小衰减和过渡带宽度的要求选择合适的窗函数类型和窗口长度N。各种类型的窗函数可达到的阻带最小衰减和过渡带宽度见表1。

表1 各种窗函数的基本参数

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

这样选定窗函数类型和长度N之后，求出单位脉冲响应

，并按照式（11）求出

。

是否满足要求，要进行演算。一般在

尾部加零使长度满足2的整数次幂，以便用FFT计算

。如果要观察细节，补零点数增多即可。必要时检查频率响应

，进行验证。

（2）窗函数法设计FIR滤波器设计实现提示

首先分析带滤波信号的频谱，根据需求，确定通带截止频率、阻带截止频率和阻带最小衰减，然后根据阻带最小衰减确定选定窗函数。

在设计之前要先写出过程，然后对照下面程序理解其实现过程。

clear

wB=1.8π; %例如矩形窗过渡带宽为1.8π/N，这里要根据选用的窗函数给出

%一方面通带截止频率和阻带截止频率可以由模拟域频率给出，转换成数字域频率

% fp=1000;  fst=2000;  Fs=10000;%模拟域频率Hz

% wp=2*fp/Fs*pi; wst=2*fst/Fs*pi; %数字域角频率弧度

% 另一方面也可以直接给出，例如

wp=0.2*pi; wst=0.3*pi; %数字域角频率

B=wst-wp; %确定过渡带宽

N=ceil(wB/B); %确定滤波器阶次

wc=(wp+wst)/2/pi;%确定滤波器3dB截止频率,归一化

hn=fir1(N-1,wc，window);% 可选窗设计低通滤波器

% 如果window缺省，默认Hamming窗，也可选择其它窗函数，具体窗函数类型可到window函数中查询。若hn=fir1(N-1,wc,boxcar(N));%选的矩形窗。hn=fir1(N-1,wc,boxcar(N));%

L=10*N;

% Hnf=fftshift(abs(fft(hn,4*N)));%幅值谱

% plot(-pi:2*pi/4/N:pi-2*pi/4/N,Hnf,'linewidth',2)；

Hnf=20*log10(abs(fft(hn,L)));%幅值谱db

b=0:2*pi/L:pi-2*pi/L;

plot(b,Hnf(1:L/2),'linewidth',2);

grid

L取不同值时，幅值谱并无明显的变化，与L相关的是频率的变化，是采样的频率

（3）滤波器设计实验

1）已知wp=0.2*pi; wst=0.3*pi

经计算得到：

Δw=wst-wp=0.1*pi; 选用矩形滤波器

     Δw=1.8*pi/N，N=18，取N=19,a=(N-1)/2, a=9

与Matlab计算结果一致。

图18    矩形窗函数幅值谱图

图19       模拟滤波器的幅值谱

2）窗函数设计的滤波器能起到滤波作用，但会有较多的噪声。

3）最终的系统函数不同，窗函数生成的系统函数中带有sa信号，该信号关于y轴对称，且信号幅度随时间上下翻转递减，导致了窗函数法生成的滤波器会有一个“尾巴”，导致滤波效果不好，模拟滤波器则没有。

2.4 语音信号滤波的滤波器设计与实现

2.4.1 加噪语音信号的滤波

1）加噪信号中有用的信号的频谱主要分布在0到4000HZ，噪声信号的频率较集中，主要有两个地方，且幅度较高；叠加之后的噪声信号频谱基本上没变换，只在两个比较集中的频率上的幅值较高。设计的滤波器应在使用后滤除掉频率在如图22中幅值最高的频率，根据实际的滤波结果选择wp=0.1pi，ws=0.2pi；δ1=3，δ2=81滤波结果最好，滤波结果如图25所示，当δ2=50时，滤波结果如图26所示。

2）根据给出的wp，ws等已知条件利用buttord函数计算N和wc；再利用butter函数计算滤波器的相关图像参数，；再利用filter函数对加噪语音信号进行滤波处理，利用fft快速傅里叶变换绘画出滤波后的频谱图。滤波器的幅度谱如图26所示。

3）设计完成后的滤波器可以满足滤波要求，虽然噪声信号仍存在，但主观感觉上和原始信号基本相同，时域信号与频域信号也基本相同。

原始信号的信噪比：SNR=   -11.9898

滤波后的信号的信噪比：SNR=   -11.9950

δ2的参数大小	SNR
原始信号	-11.9898
81	-11.9950
50	-11.9469

图20    原始语音信号时域图

图21    原始语音信号频谱图

图22    加噪语音信号频谱图

图23    加噪语音信号时域图

图24     滤波器进行滤波后的语音信号时域波形

图25   滤波处理后的语音信号的频谱图（δ2=81）

图26   滤波处理后的语音信号的频谱图(δ2=50)

图27    低通滤波器的频谱图

2.4.2 实际带噪声语音信号滤波

1） 原始带噪语音信号如图28所示，语音信号的频率在0到4KHZ之间,但主要的信号频率在2KHZ以下，2KHZ到4KHZ之间主要是高频电流产生的噪声, 根据实际的滤波结果选择wp=0.1pi，ws=0.2pi；δ1=3，δ2=50,滤波效果最好，滤波结果如图29所示，当δ2=30时滤波结果如图30所示。

图28   原始带噪语音信号频谱图

图29   滤波处理后的语音信号的频谱图

图30    滤波处理后的语音信号的频谱图

2） 根据给出的wp，ws等已知条件利用buttord函数计算N和wc；再利用butter函数计算滤波器的相关图像参数，；再利用filter函数对加噪语音信号进行滤波处理，利用fft快速傅里叶变换绘画出滤波后的频谱图。滤波器的幅度谱如图30所示。

图31     滤波器的频谱图

3）由于原始信号为真实语音信号，主观认为高频电流噪声为主要噪声，不知道电流噪声完全去除后的语音信号，无法求相对公平的SNR，但主观来说，当δ2=50 时，电流噪声已经听不到了。