滤波与功率谱估计

清华大学
《数字信号处理》期末作业
2013 年 1 月
第一题掌握去噪的方法
1.1 题目描述
MATLAB 中的数据文件noisdopp 含有噪声，该数据的抽样频率未知。

调出该数据，用你学过的滤波方法和奇异值分解的方法对其去噪。

要求：1．尽可能多地去除噪声，而又不损害原信号；
2．给出你去噪的原理与方法；给出说明去噪效果的方法或指标；
3．形成报告时应包含上述内容及必要的图形，并附上原程序。

1.2 信号特性分析
MATLAB所给noisdopp信号极其频域特征如图1.1、图1.2。

图1.1含有噪声的noisdopp信号
图1.2 noisdopp 信号频域特性
其中横坐标f 采用归一化频率，即未知抽样频率Fs 对应2（与滤波器设计时参数一致）。

信号基本特性是一个幅值和频率逐渐增加的正弦信号叠加噪声，噪声为均匀的近似白噪声，没有周期等特点。

因为噪声信号能量在全频带均匀分布，滤波器截止频率过低则信号损失大，过高则噪声抑制小，认为频谱中含有毛刺较多的部分即为信噪比较小的部分，滤除这部分可以达到较好的滤波效果。

先给定去噪效果的评定指标。

信号开始阶段频率较高（如图1.3，红圈为信号值），一周期内采样点4~5个，即信号归一化频率达到0.4~0.5（Fs=2），难以从频域将这部分信号同噪声分离，滤波后信号损失较大，故对前128点用信噪比考察其滤波效果，定义：
22()
10lg (()())k k x
k SN R y k x k =-∑∑
其中，()x k 为原nosidopp 信号，()y k 为滤波后信号。

SNR 越大表示滤除部分能力越小，可以反映滤波后信号对原信号的跟踪能力，对前128点主要考察SNR ，越大滤波器性能越好。

对128点以后的点，信号能量集中在低频部分，滤波后效果不仅仅需考察SNR ，同时考察滤波后信号的平滑程度。

定义平滑指标[1]
2
2((1)())
((1)())k k y k y k r x k x k +-=+-∑∑
r 越小表示滤波后信号相比原信号更加平滑，滤波效果更好。

图1.3 信号前100点
本文采用整体滤波、分时段滤波和奇异值分解的方法分别对信号进行去噪，并采用信噪比和平滑指数两个指标评定去噪效果。

1.3 滤波器设计
根据信号的频谱，设计IIR 滤波器和FIR 滤波器对信号去噪，用重叠相加法对信号进行卷积运算（每128点截为16段）。

根据信号的频谱，均采用低通滤波方法去噪，将滤波器截止频率设为0.2。

1.3.1 IIR 滤波器滤波
设计巴特沃斯低通滤波器，给定参数
0.17,0.23,0.1,60p s p s R dB R dB ωω====，由buttord 得到阶数7N =。

用MATLAB
函数butter 设计滤波器，并得到70点频率响应如图1.4，滤波器冲击响应如图1.5。

图1.4 巴特沃斯低通滤波器幅频特性
图1.5 巴特沃斯低通滤波器冲击响应
采用重叠相加法计算nosidopp信号于巴特沃斯低通滤波器卷积结果。

nosidopp信号每128点分为一段，共16段，得到滤波后信号有1093点。

由图1.5可以得到，设计得到的巴特沃斯滤波器存在30点的延时，取卷积后信号31点开始的1024点为滤波后信号y。

得到滤波后信号如图1.6。

与原信号比较如图1.7。

图1.6滤波后信号
图1.7滤波后信号与原信号
直观看滤波后信号去除了大部分噪声，但开始高频部分信号损失较多（如图1.7），后尾部信号频率较低，滤波后噪声含量仍比较多。

图1.8信号起始部分
前128点信噪比SNR为3.429，后896点信号SNR为12.715，全信号平滑指标为0.0642。

本题中设计的IIR滤波器的问题在于①不是线性相位，只能近似得到滤波后信号的延时，滤波后信号难以同原信号“对齐”②冲击响应无限长，有限长序列通过IIR系统得到的输出也是无限长，存在截断误差③分时段设计滤波器时计算较慢。

1.3.2 窗函数法设计FIR滤波器滤波
选择性能较优的汉宁窗设计低通滤波器，截止频率0.2，N取34，得到35点长度的滤波器。

得到滤波器幅频响应和冲击响应如图1.9，图1.10
图1.9低通滤波器幅频特性
图1.10低通滤波器冲击响应
该低通滤波器具有很好的线性相位，输出信号存在N/2=17点的延时，得到的滤波后信号如图1.11。

与原信号比较如图1.12。

图1.11滤波后信号
图1.12滤波后信号与原信号
直观感觉比巴特沃斯滤波器的滤波结果在信号尾部纹波减小，信号开始部分损失也比较大。

得到前128点信噪比SNR为6.448，后896点信号SNR为14.418，全信号平滑指标为0.0643。

由计算得到的SNR和平滑指标看，窗函数法得到的滤波器滤波效果比巴特沃斯滤波器要好。

采用矩形窗设计同样阶数的FIR滤波器，得到的滤波效果为：前128点信噪比SNR为10.303，后896点信号SNR为15.738，全信号平滑指标为0.184，如图1.13。

图1.13滤波后信号与原信号
显然，全信号采用同样的滤波器参数滤波，信号高频分量保持较好，则会保留较多的噪声。

若在频域对信号直接加窗，相当于N=∞的窗函数滤波器，恢复信号有明显的吉布斯现象，如图1.14。

图1.14滤波后信号的吉布斯现象
1.3.3 一致逼近法设计FIR滤波器滤波
选择性切比雪夫一致逼近法设计低通滤波器，截止频率0.2，N取34，通带、阻带加权均取1，得到35点长度的滤波器。

滤波器幅频响应如图1.15。

图1.15低通滤波器幅频特性
具有很好的线性相位，同样，输出信号存在N/2=17点的延时，得到的滤波后信号如图1.16。

与原信号比较如图1.17。

图1.16滤波后信号
图1.17滤波后信号与原信号
前128点信噪比SNR 为6.167，后896点信号SNR 为13.987，全信号平滑指标为0.0655，与窗函数法得到的滤波效果类似。

对全信号采用同样的滤波器参数滤波，因为没有考虑信号频率随时间的变化，不能得到更好的滤波效果。

1.4 分时段滤波
用短时傅里叶变化（128,(128),127,2nfft window hanning overlap Fs ====）估计信号的时频特性大致如图1.18。

图1.18 Nosidopp 信号的时频特性
同样将信号分为16段，每段64点，对每段信号做离散傅里叶变化求其频谱，得到幅值最大点频率为m ax f ，每两点间2/640.0313f ∆==，因为信号的包络为低频分量，需要保留，故设计m ax 3f f +∆为截止频率的低通滤波器。

采用切比雪夫一致逼近法，N 取34，通带、阻带加权均取1，输出信号存在N/2=17点的延时。

得到的滤波后信号如图1.19。

与原信号比较如图1.20。

图1.19滤波后信号
图1.20滤波后信号与原信号
图1.21滤波后信号与原信号前100点
前128点信噪比SNR为8.144，后896点信号SNR为13.734，全信号平滑指标为0.0773，得到了兼顾信号高频分量和去噪的效果，滤波效果比不分时段要好。

1.5 奇异值分解方法滤波
奇异值分解滤波的基本原理是通过矩阵分解方法，得到信号的特征值，特征值中较大的k个分量反应了信号低频成分，而较小的分量反应高频噪声，将较小的特征值置零再恢复信号，得到去噪的信号。

奇异值分解方法滤波需要考虑两个问题①一维信号构成用于奇异值分解的矩阵的方法②如何去除特征值。

本文采用教材式（9.5.18）给定的方式构造矩阵，对长度N的信号x(n)，定义矩阵X1为：
1
(0)(1)(1)
(1)(2)()
(1)()(1)
x x x M
x x x M
x L x L x N
-
⎡⎤
⎢⎥
⎢⎥=
⎢⎥⎢⎥
--
⎣⎦
X
得到的矩阵为方阵可最大程度利用其特征值，对nosidopp信号，N=1024，故选择L=513，M=512构造矩阵。

利用svd分解得到特征值，如图1.21。

同时得到矩阵U，V。

图1.22 矩阵特征值
可以看出特征值存在明显的转折点，将转折点（第30点）以后的特征值值为0。

恢复信号矩阵：
1112111
111ˆˆˆ(0)(1)(1)ˆˆˆ(1)(2)()ˆˆˆˆ(1)()(1)T x
x x
M x x x
M x
L x
L x
N -⎡⎤⎢
⎥
⎢⎥==⎢⎥⎢
⎥
--⎣⎦X U SV
此时得到的矩阵1
ˆX 将不再对称，按上式中对应位置得到滤波后信号ˆ()y x n =为：
1ˆˆ()(),0,1, (1)
k
x
i x i i N k
==-∑
得到的滤波后信号如图1.22。

与原信号比较如图1.23。

图1.23 滤波后信号
图1.24滤波后信号与原信号
前128点信噪比SNR为6.145，后896点信号SNR为13.827，全信号平滑指标为0.0708。

直观观察滤波后信号质量好与滤波器滤波后的信号。

特别对信号前100点，如图1.24，直观观察要好于同样信噪比下滤波器去噪的效果。

图1.25滤波后信号与原信号前100点
若采用将特征值中小于均值的数置零的方式滤波，得到滤波后信号如图1.26。

图1.26滤波后信号
前128点信噪比SNR为10.83，后896点信号SNR为15.357，全信号平滑指标为0.315。

对噪声去除减少，对原信号特征保持较好。

构造矩阵时，减小矩阵的行数、增加矩阵列数，或采用文献[2]中提出的构造矩阵的方法：
1
(0)(1)(1)
(1)(2)()
(1)(0)(2)
x x x M
x x x M
x N x x N
-
⎡⎤
⎢⎥
⎢⎥=
⎢⎥⎢⎥
--
⎣⎦
X
不能使得到的特征值拐点给清晰。

1.6 小结
在以信噪比SNR和输出信号平滑指标两者作为滤波效果评定标准下，采用切比雪夫一致逼近法分时段设计低通滤波器得到了最好的滤波效果，最优指标为：前128点信噪比SNR为8.144，后896点信号SNR为13.734，全信号平滑指标为0.0773。

直观观察奇异值分解去噪后滤波效果最佳。

本文去噪方式存在的不足有①给定的评定滤波效果的标准不够合理，缺乏标准信号的情况下，信噪比作为标准难以区分估计误差和噪声，给定的标准不能完
全反应直观观察出的滤波效果②缺乏对于信号所分段数L 的研究③奇异值分解方法中，对矩阵构造方式对滤波效果的研究不足。

第一题滤波器去噪部分MATLAB 程序见附录I ，奇异值分解去噪部分MATLAB 程序见附录II 。

第二题 掌握功率谱估计的方法
2.1 试验数据的产生
2.1.1 产生x(n)
第一步，产生均值为0，功率为0.12的白噪声序列12(),()u n u n 。

利用教材所给出的式(1.10.2)，可以得到给定功率的白噪声序列，即rand 函数产生的500点均
匀分布的白噪声序列，再乘以常数a ==得到12(),()u n u n 。

第二步，通过MATLAB 函数sos2tf 得到5个FIR 子系统级联组成的FIR 系统传递函数的系数firb 和fira ，firb 即为FIR 系统的冲击响应h(n)。

将第一步得到的白噪声信号12(),()u n u n 分别与h(n)卷积，即得到通过FIR 系统后的输出
12(),()v n v n ，并截取其中256
点数据得到信号12()()()v n v n jv n =+。

第三步，产生四个复正弦信号exp(2'),1,2,3,4,0,1,..,255i i a j f n i n π==，相加，再与v(n)相加，得到x(n)。

四个复正弦信号叠加得到的信号幅值和相位图如图2.1，x(n) 幅值和相位图如图2.2。

四个复正弦信号叠加信号是周期信号，周期为276。

图2.1 复正弦信号波形
图2.2 x(n)波形
2.1.2 求信号真是功率谱()j x P e ω 第一步，计算白噪声信号功率谱。

2
12()0.12(),j j u u u P e
P e ω
ω
σωπ===≤
第二步，计算FIR 系统频率响应。

()j j k
FIR k k
H e
b e
ω
ω-=
∑
第三步，计算v(n)功率谱。

由教材式(10.3.2)给出的随机信号通过线性系统理论，得到12(),()v n v n 功率谱为
2
2
112()()()
0.12
()j j j j k
j v u FIR k v k
P e
P e
H e
b e
P e
ω
ω
ω
ωω
-===∑
求v(n)功率谱，先计算自相关函数
*
1212(){()()}()(),{()()}0
v u u r m E v n v n m r m r m E u n u n m =+=++=
在计算功率谱，得到
2
1()2()0.24
j j j k
v v k k
P e
P e
b e
ω
ω
ω-==∑
第四步，求复正弦信号功率谱。

设S(n)为4个复正弦信号的和。

4
1
()exp(2')i
i i S n a
j f n π==
∑
S(n)是确定信号，且是功率信号，相求自相关函数
4
4
22()
2*21
1
(){[][]}i i i j f n j f n m j f m
s i
i i
i i r m E a e
a e
a
e
πππ+===
=
∑∑
再求功率谱
4
2
1
()()2(2),j jw m
s s i i m i P e
r m e
a f ω
πδωπωπ
∞
-=-∞
==
=-≤∑
∑
按教材式(10.4.4)，用时域傅里叶变换模平方再除以时间长度得到的功率谱为
2
1()lim {
()}21
M
j j n
s M n M
P e
E x n e
M ω
ω-→∞
=-=+∑
S(n)是周期信号，故上式化为一周期N 点内平均，ω取0~2π等分取4096点转换为DFT
2224
2
1
4
2
1
11()()(/)
(/),0,1,...,1
N
j k
j
nk
N N
s i i n i i i i P e
x n e
a N f k M N
N
N a f k N k N ππδδ-====
=
-=-=-∑
∑
∑
若信号S(n)存在截断，信号长度N ，傅里叶变换长度M ，此时S(n)功率谱为
224
21
()(/),0,1,...,1j k
M s i i i N P e
a
f k M k M M
πδ==
-=-∑
得到信号x(n)的功率谱为
2
4
2
1
()()()0.24
2(2)j j j j k
x v x k i i k i P e
P e
P e
b e
a f ω
ω
ω
ωπδωπ-==+=+-∑∑
ω取0~2π等分取4096点，离散化的功率谱
2
24
22/1
2/()
0.24
(/)j j k
x k i i k N
k
i k N
N P e
b e
a
f k M M
ω
ωωπωπδ-====+
-∑∑
取上式中N=256，M=4096，得到真实功率谱如图2.3。

图2.3 信号x 真是功率谱
2.1.3 计算给定频率点信噪比
求出复正弦信号在四个频率点处的功率。

根据教材式(3.2.32)，幅度和频率分别为',,1,2,3,4i i a f i =的复正弦信号()exp(2')i i i S n a j f n π=的功率为
2'2
11()2(2)22j si si i i i
P P e
d a f d a π
π
ω
π
π
ωπδπωωπ
π
-
-
=
=
-=⎰⎰
用离散的功率谱计算得到的结果相同，即复正弦信号在功率集中在'i f 处，幅度为2i a 。

噪声信号在'i f 处功率为220.24u σ=。

故在频率'[0.12,0.23,0.24,0.16]i f =-处的信噪比为
2
10log(/0.24)[21.79,27.78,27.78,12.22]i i SN R a ==
信噪比都比较低。

2.2 对x(n)做功率谱估计
2.2.1 周期图法估计x(n)的功率谱
周期图法估计随机信号功率谱的基本原理是根据()x n 的N 点观察值()N x n 做傅里叶变换，得到()j N X e ω，然后取其平方，再除以N ，作为对其真是功率谱的估计，周期图谱法估计的功率谱为
2
1ˆ()()P E R N
P X N ωω=
令ω在0~2π等分取4096点，得到离散化的功率谱
2
1ˆ()()PER N
P k X k N
=
对256点()x n 做4096点DFT ，得到估计的曲线如图2.4。

图2.4 周期图法估计x 功率谱
得到比较好的功率谱估计值，256点实际分辨率为2/2560.0078<0.01B ==，故四个频率点的谱线均可分开，能量较小'0.16f =-处也有所反应。

2.2.2 Welch 方法改进周期图 W 点汉明窗表达式为
2()0.540.46cos(
),0,1,...,1n w n n W W
π=-=-
每长度M=64分一段，得到Bartlett 法估计的每段功率谱表达式为
2
1
1ˆ()()M i i j n
PER
N n P x n w n e
M U
ω--==∑
其中
1
2
1()M n U w
n M
-==
∑
平均后的功率谱作为最后估计的功率谱
1
1ˆ()L
i PER
PER
i P P L
ω==∑
令ω在0~2π等分取4096点，，用DFT 代替DTFT ，可得到离散化的功率谱估计。

取M=64，每段重叠32点，分7段，用上述Welch 法得到的功率谱估计曲线如图2.5。

图2.5 Welch 法估计x 功率谱
功率谱较周期谱法更光滑，但是分辨率下降，不能区分0.23、0.24两处频率值的信号。

64点汉明窗的主瓣宽度为04/0.06250.01B N ==>，故不能分辨两个频率值，若到达0.01的分辨率，至少需400点的窗函数，超过信号长度，等同于周期谱法，故改进的周期谱法不能达到更好的功率谱估计效果。

2.2.3 自相关法估计功率谱
自相关法估计功率谱的基本原理是，根据信号的N 点观察值()N x n ，计算自相关函数的估计值ˆ()r m ，再对ˆ()r m 做傅里叶变化得到估计的功率谱。

复信号的自相关函数ˆ()r m 计算方法为
1*
1ˆ()()()N m
N N n r
m x n x n m N
--==+∑
再得到功率谱估计为
0ˆˆˆ()2R e[()](0)M
j m BT
m P r m e r ωω-==-∑ 根据上述公式，取M=63，令ω在0~2π等分取4096点，用DFT 代替DTFT ，得到的离散的功率谱计算方式为
20ˆˆˆ()2R e[()](0)M
j
m k
N BT
m P k r m e r
π-==-∑ 得到估计的功率谱曲线如图2.6（已舍去负值）。

图2.6 自相关法估计x 功率谱
谱线比周期谱法光滑，但不能分辨0.23、0.24频率处的两个谱线，减小M 或使用窗函数可使谱线光滑，但将更加降低频率分辨率。

2.2.4 AR 模型估计功率谱的自相关法
自相关法求解AR 模型的基本原理是用信号的N 点观察值()N x n 计算得到的自相关函数估计值ˆ()r m 代替Yule-Walker 方程中的()r m ，得到AR 模型参数p a 。

采用Levinson-Durbin 递推算法计算。

初始条件
0ˆ(0)(0)x x r r
ρ== 递推关系为
1
111
112
1ˆˆ[()()()]/()()()
[1]
m m m x x m k m m m m m m m k a k r
m k r m a k a k k a m k k ρρρ---=---=--+=+-=-∑
用MATLAB 函数levinson 可以求解上述方程参数，得到p 阶AR 模型参数。

自相关求解AR 模型的最小预测误差随p 阶数增加而减小，p 为6时，最小预测误差77.15，p 为15时最小预测误差为43.3，取p 为15，得到AR 模型系数
15ˆ(),0,1,2,..,a
k k p =，和最小预测误差15ˆρ
，从而得到估计的功率谱 152
15
ˆˆ()ˆ()AR
p
j k
k P a
k e
ωρ
ω-==∑
N=4096离散得到
215221
150
ˆˆ()ˆ()j l
N AR
N j
lk
N
k P e a
k e ππρ
--==
∑
得到估计的功率谱曲线如图2.7。

图2.7 AR 模型自相关算法估计x 功率谱，p=15
基本估计出四处频率点的复正弦信号，但不能分辨0.23和0.24两处频率的谱线，信噪比较低的-0.16频率点处的信号估计不佳，若降低阶数p ，则根本估计不出该处的谱线，p 取10时得到的功率谱如图2.8。

图2.8 AR 模型自相关算法估计x 功率谱，p=10
2.2.5 AR 模型估计功率谱的Burg 算法 根据教所给的Burg 算法计算。

设定初始条件00()()f b e n x n e ==，由
1
*
111
1
2
2
112()(1)
ˆ,1,2,...,()(1)
N f b m m n m
m
N N f
b
m m n m
n m
e n e n k m p e n e n ---=----==--==+
-∑∑
∑
得到1ˆk 。

由ˆ(0)x r 得到m=1时的参数21111
ˆˆˆ(1),(1)(0)x a k k r ρ==-。

再由 11*
1100()()(1)()()()()()(),1,2,...,f f b
m m m m b
b
f
m m m m f
f
e n e n k e n e n e n k e n e n e n x n m p
----=+-=+===
求得11(),()f b e n e n ，由ˆm k 表达式求得2ˆk 。

由
*112
1ˆˆˆˆ()()(),1,2,...,1ˆˆ()ˆˆˆ(1)m m m m m m m m
m a k a k k a n k k m a m k k ρ
ρ---=+-=-==-
求得m=2时参数222ˆˆˆ(1),(2),a
a ρ。

重复上述步骤，可得到m=p 时参数。

根据上述算法，在MATLAB 中编写程序实现上述算法。

得到p=18时的参数
10ˆ(),0,1,...,18,fb
a
k k ρ=，并利用MATLAB 函数arburg 验证了计算结果的正确性。

得到估计的功率谱如图2.9。

图2.9 AR 模型Burg 算法估计x 功率谱，p=18
比自相关法有更好的频谱，可以分辨0.23和0.24频率处的谱线，-0.16处的信号也更加明显。

但Burg 算法在p 较小时（如10），得到的功率谱不仅不能分辨
较近的频谱和幅度小的信号，且谱线的相对幅度与信号幅度对应发生变化，如图2.10。

图2.10 AR 模型Burg 算法估计x 功率谱，p=10
2.2.6 AR 模型估计功率谱的Burg 算法 根据信息论准则
()ln()2k AIC k N k
ρ=+
利用自相关算法，k 由1开始增加，通过levinson 函数求解最小估计误差k ρ，计算()AIC k ，得到最小值10k =，即为最合适的阶次p ，此时(10)=1017.8AIC ，最小估计误差49.38。

p 为10时利用自相关算法的AR 模型的功率谱估计结果如前2.8图。

2.3 小结
对本题中信号的功率谱估计不足在于①产生的实验信号中的白噪声不是理想白噪声，不同次实验产生的12(),()u n u n 功率不稳定，两者的独立性也不佳②缺乏对功率谱估计效果的定量指标。

2.1~2.2.4，2.2.6节相关MATLAB 程序见附录III ，2.2.5节相关MA TLAB 程序见附录VI 。

参考文献
[1] 陈强, 黄声享, 王韦. 小波去噪效果评价的另一指标[J]. 测绘信息与工程, 2008(5):13-14.
[2] 周江, 杨浩. 分段滤波在心率变异性功率谱分析中的应用[J]. 重庆大学学报(自然科学版),
2001(4):95-97.
[3] 海兰萍, 姜占才, 李振起. 一种改进的奇异值分解语音降噪算法[J]. 科技信息,
2012(6):156-158.
[4] 迟华山, 王红星, 赵培洪, 等. 基于S变换的线性调频信号时频滤波[J]. 无线电通信技术,
2012(1):21-24.
[5] DELL. 常用图像去噪比较与分析[Z]. 20101.
[6] 奇异值分解降噪的改进方法-张磊彭伟才原春晖刘彦[Z].
[7] 胡广书. 数字信号处理[M]. 第二版. 北京: 清华大学出版社, 2008.
[8] 自适应线性调频信号时变滤波器在分数傅里叶域[Z].
[9] 姚文俊. 自相关法和Burg法在A R模型功率谱估计中的仿真研究[J]. 计算机与数字工程,
2007(10):32-34.
[10] 白噪声背景下正弦信号频率和功率估计的有效算法[Z].
[11] 李军, 王凯. 三种现代功率谱估计方法性能研究[J]. 科技信息, 2010(30):554-555.
附录I
% 期末大作业第一题
clear all;clc;close all
load noisdopp%读入信号
%信号特性
x=noisdopp;nx=length(x);
plot(x);axis tight;xlabel('n');ylabel('noisdopp');grid on
plot(1:100,x(1:100),1:100,x(1:100),'ro');axis tight;xlabel('n');ylabel('noisdopp');grid on
x_fft=fft(x);figure
subplot(2,1,1);plot([1:nx/2]*2./nx,abs(x_fft(1:nx/2)));axis tight;xlabel('f');ylabel('|X(k)|');
subplot(2,1,2);plot([1:nx/2]*2./nx,unwrap(angle(x_fft(1:nx/2))));axis tight;xlabel('f');ylabel('phase X(k)'); % 分段处理信号；每M点分一段
M=128/2;
% 设计滤波器
%FIR，窗函数法阶数为L，有L+1点，有L/2点的延时
L=34;fc=0.2; %Fs=2
b1=fir1(L,fc,hanning(L+1));figure;freqz(b1);figure;stem(b1);xlabel('n');ylabel('h(n)');grid on
%FIR，切比学夫一致逼近法
F=[0,fc-0.03,fc+0.03,1]; A=[1,1,0,0]; W=[1,1];b2=remez(L,F,A,W);freqz(b2)
%将H(1)置为一
% b2_mo=max(abs(fft(b2)));% b2=remez(L,F,A,W)./b2_mo;
%IIR，巴特沃斯滤波器
[Niir,wniir]=buttord(fc-0.03,fc+0.03,0.1,60);[bb,ab]=butter(Niir,fc);
%[bb,ab]=bilinear(bbs,abs,2);
L_b3=2*L+1;b3=impz(bb,ab,L_b3);freqz(b3);figure;stem(b3);xlabel('n');ylabel('h(n)');grid on
%FIR，窗函数法,吉布斯现象
nb51=1:L/2-1;nb52=L/2+1:L+1;wc_b5=fc*2*pi;
b5=[sin((nb51-L/2)*wc_b5)./(pi*(nb51-L/2)),wc_b5./pi,sin((nb52-L/2)*wc_b5)./(pi*(nb52-L/2))]; % 用于重建信号
y1=zeros(1,nx+length(b1)-1);y2=zeros(1,nx+length(b2)-1);y3=zeros(1,nx+length(b3)-1);
y4=zeros(1,nx+length(b1)-1);y5=zeros(1,nx+length(b5)-1);ylw=zeros(1,nx);
for n=1:(nx-M)/(M-n_over)+1
xw=x((n-1)*(M-n_over)+1:(n-1)*(M-n_over)+M);
%直接频域加窗
xw_fft=fft(xw); xl=[xw_fft(1:5),zeros(1,length(xw)-10),xw_fft(length(xw)-4:length(xw))];
yl=ifft(xl,'symmetric');
ylw((n-1)*(M-n_over)+1:(n-1)*(M-n_over)+M)=ylw((n-1)*(M-n_over)+1:(n-1)*(M-n_over)+M)+yl;
%分段调整滤波频率,取频率最大值左右一定带宽的信号
[max_fft,n_max_fft]=max(xw_fft(1:M/2));
%通带中心频率
f0=n_max_fft./M*2;
%设置带宽
fw=3;df=1./M*2;
if n_max_fft-fw<=0
%即为低通滤波器，FIR，切比学夫一致逼近法
F=[0,f0,f0+fw*df,1]; b4=remez(L,F,A,W);
else
%F=[0,f0-fw*df,f0-fw*df+df,f0+fw*df-df,f0+fw*df,1]; %即为带通滤波器
F=[0,f0,f0+fw*df,1];b4=remez(L,F,A,W);
end
%信号还原
yw5=conv(xw,b5);
y5((n-1)*(M-n_over)+1:(n-1)*(M-n_over)+M+length(b5)-1)=y5((n-1)*(M-n_over)+1:(n-1)*(M-n_over)+M+ length(b5)-1)+yw5;
yw4=conv(xw,b4);
y4((n-1)*(M-n_over)+1:(n-1)*(M-n_over)+M+length(b4)-1)=y4((n-1)*(M-n_over)+1:(n-1)*(M-n_over)+M+ length(b4)-1)+yw4;
yw3=conv(xw,b3);
y3((n-1)*(M-n_over)+1:(n-1)*(M-n_over)+M+length(b3)-1)=y3((n-1)*(M-n_over)+1:(n-1)*(M-n_over)+M+ length(b3)-1)+yw3.';
yw2=conv(xw,b2);
y2((n-1)*(M-n_over)+1:(n-1)*(M-n_over)+M+length(b2)-1)=y2((n-1)*(M-n_over)+1:(n-1)*(M-n_over)+M+ length(b2)-1)+yw2;
yw1=conv(xw,b1); %信号幅值不对
y1((n-1)*(M-n_over)+1:(n-1)*(M-n_over)+M+length(b1)-1)=y1((n-1)*(M-n_over)+1:(n-1)*(M-n_over)+M+l ength(b1)-1)+yw1;
end
figure
[my1,ny_y1]=max(b1);y1_yx=y1(ny_y1:length(y1));plot(y1_yx);axis tight;xlabel('n');ylabel('y');grid on figure;plot(1:nx,x,1:nx,y1_yx(1:nx),'r-');axis tight;xlabel('n');ylabel('x,y');grid on
[my2,ny_y2]=max(b2);y2_yx=y2(ny_y2:length(y2));
plot(y2_yx);axis tight;xlabel('n');ylabel('y');grid on
figure;plot(1:nx,x,1:nx,y2_yx(1:nx),'r-');axis tight;xlabel('n');ylabel('x,y');grid on
%figure;%plot(ylw) % 直接频域加窗，等于时域用sinc函数截断，效果不好，吉布斯现象
[my3,ny_y3]=max(b3);y3_yx=y3(ny_y3:length(y3));plot(y3_yx);axis tight;xlabel('n');ylabel('y');grid on plot(1:nx,x,1:nx,y3_yx(1:nx),'r-');axis tight;xlabel('n');ylabel('x,y');grid on
[my4,ny_y4]=max(b4);y4_yx=y4(ny_y4:length(y4));plot(y4_yx);axis tight;xlabel('n');ylabel('y');grid on figure;plot(1:nx,x,1:nx,y4_yx(1:nx),'r-');axis tight;xlabel('n');ylabel('x,y');grid on
[my5,ny_y5]=max(b5);y5_yx=y5(ny_y5:length(y5));% plot(y5_yx);axis tight;xlabel('n');ylabel('y');grid on % plot(1:nx,x,1:nx,y5_yx(1:nx),'r-');axis tight;xlabel('n');ylabel('x,y');grid on
%计算峰值信噪比
% y_MSN=y3_yx(129:nx);% y_MSN128=y3_yx(1:128);% y_MSN=y1_yx(129:nx);
% y_MSN128=y1_yx(1:128);% y_MSN=y5_yx(129:nx);% y_MSN128=y5_yx(1:128);
% y_MSN=y2_yx(129:nx);% y_MSN128=y2_yx(1:128);
y_MSN=y4_yx(129:nx);y_MSN128=y4_yx(1:128);
x_MSN128=x(1:128);x_MSN=x(129:nx);
MSN=(y_MSN-x_MSN)*(y_MSN-x_MSN)';
PSNR=10*log10(x_MSN*x_MSN'/MSN);
MSN128=(y_MSN128-x_MSN128)*(y_MSN128-x_MSN128)';
PSNR128=10*log10(x_MSN128*x_MSN128'/MSN128);
%平滑尺度
% y_r1=y3_yx(2:nx);y_r2=y3_yx(1:nx-1);% y_r1=y1_yx(2:nx);y_r2=y1_yx(1:nx-1);
% y_r1=y5_yx(2:nx);y_r2=y5_yx(1:nx-1);% y_r1=y2_yx(2:nx);y_r2=y2_yx(1:nx-1);
y_r1=y4_yx(2:nx);y_r2=y4_yx(1:nx-1);
x_r1=x(2:nx);x_r2=x(1:nx-1);r_phx=(x_r1-x_r2)*(x_r1-x_r2)';
r_phy=(y_r1-y_r2)*(y_r1-y_r2)';r_ph=r_phy/r_phx;
% 用短时傅里叶变化估计信号时频特性
% spectrogram(x,128,127,128,2);
%计算自相关函数
% rx=xcorr(x,x);subplot(3,1,1);plot(rx);xlabel('n');ylabel('r')
% ry=xcorr(y,y);subplot(3,1,2);plot(ry);
% ry1=xcorr(y1,y1);subplot(3,1,3);plot(ry1);
%对不足一周期正弦信号FFT分析效果
% nt=0:1/300:127/300;
% yt=sin(2*pi*nt);
% Eb=(b4*b4');
% ytl=conv(yt,b4);
% plot([yt',ytl(1:length(yt))'])
%plot(abs(fft(yt)))
附录II
% 奇异值分解方法降噪
clear all;clc;close all;load noisdopp%读入信号
%信号特性
x=noisdopp;%plot(x);
nx=length(x);
% LxL=nx方法
M=nx/2;X1=zeros(M+1,M);
for i=1:M+1
for j=1:M
X1(i,j)=x(i+j-1);
end
end
[U,S,V]=svd(X1);plot(abs(diag(S)));axis tight;xlabel('k');grid on
%滤波
S1=zeros(i,j);
% S1(1:30,1:30)=S(1:30,1:30); %衰减很明显
s_u=sum(diag(S))./length(diag(S));s_n=find(diag(S)-s_u>0);s_nq=max(s_n);S1(1:s_nq,1:s_nq)=S(1:s_nq,1:s_ nq); %均值下去除
%重建,奇异值分解降噪的改进方法
A=U*S1*V.';y=0*x;
for i=1:nx
ytemp=0;
if i<=M
for j=1:i
ytemp=ytemp+A(j,i-j+1);
end
y(i)=ytemp/i;
else
for j=1:nx-i+1
ytemp=ytemp+A(M+1-j+1,i-M+j-1);
end
y(i)=ytemp/(nx-i+1);
end
end
plot(y);axis tight;xlabel('n');ylabel('y');grid on;figure
plot(1:nx,x,1:nx,y,'r-');axis tight;xlabel('n');ylabel('x,y');grid on
%计算峰值信噪比
y_MSN=y(129:nx);y_MSN128=y(1:128);x_MSN128=x(1:128);x_MSN=x(129:nx);
MSN=(y_MSN-x_MSN)*(y_MSN-x_MSN)';PSNR=10*log10(x_MSN*x_MSN'/MSN);
MSN128=(y_MSN128-x_MSN128)*(y_MSN128-x_MSN128)';
PSNR128=10*log10(x_MSN128*x_MSN128'/MSN128);
%平滑尺度
y_r1=y(2:nx);y_r2=y(1:nx-1);x_r1=x(2:nx);x_r2=x(1:nx-1);r_phx=(x_r1-x_r2)*(x_r1-x_r2)';
r_phy=(y_r1-y_r2)*(y_r1-y_r2)';r_ph=r_phy/r_phx;
附录III
% 第二题
clear all;clc;close all
N=500; %信号长度
%功率为0.12的白噪声产生方式(1.10.2)
P=0.12;a=sqrt(12*P);
%产生噪声序列
u1=rand(1,N)*a;u2=rand(1,N)*a;
%5个FIR子系统
B1=[1,1.98,0.98];B2=[1,-1.98,0.98];B3=[1,-1.8418,0.98];B4=[1,-1.5,0.98];
B5=[1,-1.2727,0.98];B=[B1;B2;B3;B4;B5];sos=[B,ones(5,1),zeros(5,2)];
%得到级联系统传递函数
[firb,fira]=sos2tf(sos);v1=conv(u1,firb);v2=conv(u2,firb);nx=256;
v=v1(100:100+nx-1)+1i*v2(100:100+nx-1);
%产生复正弦信号
exp_a=[6,12,12,2];exp_f=[0.12,0.23,0.24,-0.16];
for i=1:length(exp_a)
exp_n=0:nx-1;
exp_x(i,1:nx)=exp_a(i)*exp(1i*2*pi*exp_f(i)*exp_n);
end
exp_s=sum(exp_x);x=exp_s+v;
save('x_sv.mat','x') %保存用于burg算法
subplot(2,1,1);plot(abs(exp_s));axis tight;xlabel('n');ylabel('|exp|');grid on
subplot(2,1,2);plot(angle(exp_s));axis tight;xlabel('n');ylabel('phase exp/rad');grid on %画x波形图
figure
subplot(2,1,1);plot(abs(x));axis tight;xlabel('n');ylabel('|x|');grid on
subplot(2,1,2);plot(angle(x));axis tight;xlabel('n');ylabel('phase x/rad');grid on
M=4096;nw=-M/2:M/2-1;fw=nw./M;w=nw*(2*pi./M);bt=zeros(1,length(w)); for k=0:10
bt=firb(k+1).*exp(-1i*w*k)+bt;
end
Pv=abs(bt).^2*2*P;%Ps=1./length(w)*(abs(fft(exp_s,length(w))).^2);
Ps=zeros(1,length(w));
for i=1:4
expndt=find(2*pi*exp_f(i)-w<0); expd=min(expndt);
Ps(expd)=length(x).^2./length(w)*exp_a(i).^2;
end
% plot(Ps)
Px=Pv+Ps;figure;plot(fw,Px);grid on;axis tight;xlabel('f');ylabel('Px');
%计算信噪比
SNR=10*log10(exp_a.^2./2/P);
%周期图法功率谱估计
%求信号FFT
Xper=fftshift(fft(x,M));Pper=abs(Xper).^2./M;figure
plot(fw,Pper);grid on;axis tight;xlabel('f');ylabel('Pper');
%用Welch法估计
mw=64;nxw=length(x);lw=(nxw-mw/2)/(mw/2);
%得到汉明窗
wnw=0:mw-1;wn_h=0.54-0.46*cos(2*pi*wnw/mw);Uw=1./mw*(wn_h*wn_h');
%将x分段
Pper_w=zeros(1,M);
for i=1:lw
xni=(x((i-1)*mw/2+1:(i+1)*mw/2)).*wn_h; xni_fft=fft(xni,M);
Pwt=1./(M*Uw)*(abs(xni_fft).^2); Pper_w=Pper_w+Pwt;
end
Pper_w=1./lw*fftshift(Pper_w)*16;figure;plot(fw,Pper_w);grid on;axis tight;xlabel('f');ylabel('Pper Welch'); %自相关法估计
m_bt=63;n_bt=length(x);
%计算自相关函数
r_bt=zeros(1,m_bt);
for i=0:m_bt-1
r_bt(i+1)=1./n_bt*sum(x(1:n_bt-i)'.'.*x(1+i:n_bt));
end
Pbt=(2*real(fftshift(fft(r_bt,M)))-r_bt(1))./2/pi;figure;
plot(fw,Pbt);grid on;axis([-0.5 0.5 0 max(Pbt)]);xlabel('f');ylabel('Pper 自相关法');
%AR模型估计%自相关法
p=10;[a_ar,E_ar,k_ar]=levinson(r_bt(1:p+1),p);
%求功率谱
a_arf=fftshift(fft(a_ar,M));Par=E_ar./(abs(a_arf).^2)./32;figure
plot(fw,Par);grid on;axis tight;xlabel('f');ylabel('Par');
%AIC法确定最佳阶次
k=2;[a_ar,E_ar,k_ar]=levinson(r_bt(1:1+1),1);aic(1)=n_bt*log(E_ar)+2;
[a_ar,E_ar,k_ar]=levinson(r_bt(1:2+1),2);aic(2)=n_bt*log(E_ar)+4;
while aic(k)<aic(k-1) && k<=30
k=k+1; [a_ar,E_ar,k_ar]=levinson(r_bt(1:k+1),k); aic(k)=n_bt*log(E_ar)+2*k;
end
%pburg(x,10,M)
%P=1./length(w)*(abs(fft([v,v,v],length(w))).^2);
% figure% plot(Pv);%Pv1=abs(firb(2).*(exp(-1i*w)+1i*exp(-1i*w))+firb(1)*(1+1i)).^2;
% figure% plot(P);
附录VI
% burg算法
clc;close all;clear all
%读入信号，同前算法信号
load('x_sv.mat');p=10;N=length(x);m=1;a=zeros(p,p);ef=zeros(p,N);
%设定初始条件
ef0=x;eb0=x;
k(m)=(-2*sum(ef0(2:N).*eb0(1:N-1)'.'))./(ef0(2:N)*ef0(2:N)'+eb0(1:N-1)*eb0(1:N-1)'); rox=1./N*(x*x');a(m,1)=k(m);ro(m)=(1-k(m)*k(m)')*rox;
ef=ef0(2:N)+k(m)*eb0(1:N-1);eb=eb0(1:N-1)+k(m)'*ef0(2:N);
while m<p
m=m+1;efm=ef(2:length(ef));ebm=eb(1:length(eb)-1);
k(m)=-2.*efm*ebm'./(efm*efm'+ebm*ebm');
ef=efm+k(m)*ebm;eb=ebm+k(m)'*efm;
for ak=1:m-1
a(m,ak)=a(m-1,ak)+k(m)*a(m-1,m-ak)';
end
a(m,m)=k(m);ro(m)=(1-k(m)*k(m)')*ro(m-1);
end
%[ab,Eb,kb]=arburg(x,p) %验证计算结果
%求功率谱
ropf=1./(N-p)*(ef*ef');ropb=1./(N-p)*(eb*eb');rofb=1/2*(ropf+ropb);
M=4096;nw=-M/2:M/2-1;fw=nw./M;w=nw*(2*pi./M);
a_arf=fftshift(fft([1,a(p,:)],M));Par=rofb./(abs(a_arf).^2)./128;
figure;plot(fw,Par);grid on;axis tight;xlabel('f');ylabel('Par');。