滤波与功率谱估计
功率谱计算[解说]
![功率谱计算[解说]](https://img.taocdn.com/s3/m/3cec2a18f02d2af90242a8956bec0975f465a478.png)
功率谱计算功率谱估计在现代信号处理中是一个很重要的课题，涉及的问题很多。
在这里，结合matlab，我做一个粗略介绍。
功率谱估计可以分为经典谱估计方法与现代谱估计方法。
经典谱估计中最简单的就是周期图法，又分为直接法与间接法。
直接法先取N点数据的傅里叶变换（即频谱），然后取频谱与其共轭的乘积，就得到功率谱的估计；间接法先计算N点样本数据的自相关函数，然后取自相关函数的傅里叶变换，即得到功率谱的估计.都可以编程实现，很简单。
在matlab中，周期图法可以用函数periodogram实现。
但是周期图法估计出的功率谱不够精细，分辨率比较低。
因此需要对周期图法进行修正，可以将信号序列x(n)分为n个不相重叠的小段，分别用周期图法进行谱估计，然后将这n段数据估计的结果的平均值作为整段数据功率谱估计的结果。
还可以将信号序列x(n)重叠分段，分别计算功率谱，再计算平均值作为整段数据的功率谱估计。
这2种称为分段平均周期图法，一般后者比前者效果好。
加窗平均周期图法是对分段平均周期图法的改进，即在数据分段后，对每段数据加一个非矩形窗进行预处理，然后在按分段平均周期图法估计功率谱。
相对于分段平均周期图法，加窗平均周期图法可以减小频率泄漏，增加频峰的宽度。
welch法就是利用改进的平均周期图法估计估计随机信号的功率谱，它采用信号分段重叠，加窗，FFT等技术来计算功率谱。
与周期图法比较，welch法可以改善估计谱曲线的光滑性，大大提高谱估计的分辨率。
matlab中，welch法用函数psd实现。
调用格式如下：[Pxx,F] = PSD(X,NFFT,Fs,WINDOW,NOVERLAP)X：输入样本数据NFFT：FFT点数Fs:采样率WINDOW：窗类型NOVERLAP，重叠长度现代谱估计主要针对经典谱估计分辨率低和方差性不好提出的，可以极大的提高估计的分辨率和平滑性。
可以分为参数模型谱估计和非参数模型谱估计。
参数模型谱估计有AR模型，MA模型，ARMA模型等；非参数模型谱估计有最小方差法和MUSIC法等。
matlab中 功率谱估计的函数
在matlab中，功率谱估计是信号处理和频谱分析中常用的一种方法。
通过对信号的频谱特性进行估计，可以有效地分析信号的功率分布情况，从而为信号处理和系统设计提供重要的参考信息。
在matlab中，提供了多种功率谱估计的函数，以下将对其中几种常用的函数进行介绍和分析。
1. periodogram函数periodogram函数是matlab中用于估计信号功率谱密度的函数之一。
它基于傅里叶变换将离散时间信号转换成频域信号，然后计算频域信号的功率谱密度。
其调用格式为：[Pxx, F] = periodogram(x,window,nfft,fs)其中，x为输入的离散时间信号，window为窗函数，nfft为离散傅里叶变换的点数，fs为信号的采样频率。
periodogram函数返回的Pxx 为功率谱密度估计值，F为对应的频率。
2. pwelch函数pwelch函数也是用于估计功率谱密度的函数，它采用了Welch方法，通过对信号进行分段处理，然后对各段信号进行傅里叶变换，并对各段功率谱密度进行平均。
其调用格式为：[Pxx, F] = pwelch(x,window,noverlap,nfft,fs)其中，x为输入的离散时间信号，window为窗函数，noverlap为相邻分段的重叠点数，nfft为离散傅里叶变换的点数，fs为信号的采样频率。
pwelch函数返回的Pxx为功率谱密度估计值，F为对应的频率。
3. cpsd函数cpsd函数用于估计信号的交叉功率谱密度，即两个信号之间的频谱特性。
其调用格式为：[Pxy, F] = cpsd(x,y,window,noverlap,nfft,fs)其中，x和y为输入的两个离散时间信号，window为窗函数，noverlap为相邻分段的重叠点数，nfft为离散傅里叶变换的点数，fs为信号的采样频率。
cpsd函数返回的Pxy为交叉功率谱密度估计值，F为对应的频率。
4. mscohere函数mscohere函数用于估计信号的相干函数，即两个信号之间的相关性。
第6讲 谱估计2.传统法
作长度为L>2M－1的序列：
RN (m) w(m) Qx (m) = 0 R (m L) w(m L) N 0 ≤ m ≤ M 1 M ≤ m ≤ L M +1 L M +1 < m ≤ L 1
求功率谱估计 FFT[Qx (m)] = S x (k 2π / L) k = 0,1, , L 1 2．估计质量 ． 当M<<N-1，BT法不等于周期图法，而是对周期图 的平滑，以下讨论BT谱或者周期图平滑谱的性能。
加窗后的功率谱估计值为：
S x (ω ) =
w(m) W (e )
jω
M 1
m = ( M 1)
w ( m ) R N ( S x (ω ) = ，有 RN (m) S (e jω ) = I N (ω ) 2π
∧
∫π
π
I N (θ )W (e j (ω θ ) )dθ
4.7.1 间接法
间接法是由随机信号N个观察值 x(0),, x( N 1) ，估计 出自相关函数 RN (m) ，然后求其傅里叶变换作为功 率谱的估计。
S (e ) =
∧ jω m= M
RN ( m )e jωm ∑
M
M ≤ N 1
4.7.2 直接法
设有限长实序列
x( n) (0 ≤ n ≤ N 1) 的傅里叶变换为：
ω
∧
∧
线性系统幅度频率响应的平方是 S x (e jω )，因此
sin[ωN ]2 2 此时谱估计方差为： [S (ω)] = S x2 (e jω )1 + Var N sin ω
∧
S (ω ) ≈ S x (e ) S ω (ω )
∧
jω
功率谱和功率谱密度计算公式
功率谱和功率谱密度计算公式
功率谱（Power Spectrum）
是描述随机信号或时间序列在不同频率下功率分布情况的工具。
对于离散信号，功率谱的计算通常涉及到傅里叶变换（Fourier Transform）或者更一般的傅里叶分析方法。
假设有一个离散信号(x(n))（其中(n)表示时间或样本序号），其功率谱(P(f))可以通过以下步骤计算：
傅里叶变换：首先，对信号(x(n))进行傅里叶变换，得到其频谱(X(f))：
(X(f) = \sum_{n=-\infty}^{\infty} x(n) e^{-j2\pi fn})
计算功率谱：然后，计算频谱的模的平方，即得到功率谱(P(f))：
(P(f) = |X(f)|^2)
功率谱密度（Power Spectral Density, PSD）
是单位频率范围内的平均功率，通常用于描述连续信号的功率分布。
对于连续信号(x(t))（其中(t)表示时间），其功率谱密度(S_{xx}(f))可以通过自相关函数和傅里叶变换得到：
自相关函数：首先，计算信号(x(t))的自相关函数(R_{xx}(\tau))：
(R{xx}(\tau) = \int{-\infty}^{\infty} x(t) x(t+\tau) dt)
傅里叶变换：然后，对自相关函数(R{xx}(\tau))进行傅里叶变换，得到功率谱密度(S{xx}(f))：(S{xx}(f) = \int{-\infty}^{\infty} R_{xx}(\tau) e^{-j2\pi f\tau} d\tau)。
自-信号检测与估计试题——答案(不完整版)
一、概念:1. 匹配滤波器。
概念：所谓匹配滤波器是指输出判决时刻信噪比最大的最佳线性滤波器。
应用:在数字信号检测和雷达信号的检测中具有特别重要的意义。
在输出信噪比最大准则下设计一个线性滤波器是具有实际意义的。
2. 卡尔曼滤波工作原理及其基本公式(百度百科)首先,我们先要引入一个离散控制过程的系统。
该系统可用一个线性随机微分方程（Linear StoｃhastiｃDifferｅｎce equatiｏn）来描述：X(k）=A X(k－1)+BＵ(k)+W(ｋ)再加上系统的测量值:Z(k)＝ＨＸ(k)+V(k)上两式子中，Ｘ(k)是k时刻的系统状态,U(k）是k时刻对系统的控制量。
A和B是系统参数，对于多模型系统,他们为矩阵。
Z(k)是k时刻的测量值,H是测量系统的参数，对于多测量系统,H为矩阵。
Ｗ(k)和V（k)分别表示过程和测量的噪声。
他们被假设成高斯白噪声(Whｉte ＧausｓｉaｎNoiｓe），他们的coｖaｒiancｅ分别是Q，R（这里我们假设他们不随系统状态变化而变化）。
对于满足上面的条件(线性随机微分系统,过程和测量都是高斯白噪声),卡尔曼滤波器是最优的信息处理器。
下面我们来用他们结合他们的covariａncｅs来估算系统的最优化输出（类似上一节那个温度的例子）。
首先我们要利用系统的过程模型，来预测下一状态的系统。
假设现在的系统状态是k,根据系统的模型,可以基于系统的上一状态而预测出现在状态：Ｘ(k|ｋ－1）=A X(k-１|k－1）+B U(k) ……….．(1)式(１)中，X（k|k-1)是利用上一状态预测的结果，X（k-１|k－１)是上一状态最优的结果,Ｕ(k)为现在状态的控制量，如果没有控制量,它可以为0。
到现在为止，我们的系统结果已经更新了,可是,对应于X(ｋ｜k-1)的covａｒiancｅ还没更新。
我们用P表示covａｒianｃe:P(k|ｋ-１)=A P(k-1｜ｋ－１)Ａ’+Q (2)式(2)中,P（k｜ｋ－1）是X（k|k-1)对应的covaｒｉaｎce,P(ｋ-１|k-1）是Ｘ(k-1|k-1)对应的cｏvariaｎce,A’表示A的转置矩阵,Q是系统过程的cｏvａriance。
经典功率谱估计
雷达和声呐系统
目标检测
在雷达和声呐系统中，经典功率谱估计常被用于目标检测。通过对接收到的信号进行功率 谱分析，可以判断是否存在目标以及目标的位置和速度等信息。
距离和速度测量
在雷达和声呐系统中，经典功率谱估计还可以用于距离和速度测量。通过对接收到的信号 进行功率谱分析，可以估计出目标与系统之间的距离和相对速度。
信号分类
在雷达和声呐系统中，经典功率谱估计还可以用于信号分类。通过对接收到的信号进行功 率谱分析，可以判断目标的类型，例如区分飞机、船舶或车辆等不同类型目标。
05 经典功率谱估计的改进方 法
基于小波变换的功率谱估计
1
小波变换能够将信号分解成不同频率和时间尺度 的分量，从而更好地揭示信号的内在结构和特征。
然而，这些方法通常需要较长 的数据长度和较为复杂的计算 过程，对于短数据和实时处理 的应用场景具有一定的局限性 。
研究展望
01
随着信号处理技术的发展，经典功率谱估计方法仍有进一步优化的空 间。
02
针对短数据和实时处理的应用场景，研究更为快速、准确的功率谱估 计方法具有重要的实际意义。
03
结合机器学习和人工智能技术，探索基于数据驱动的功率谱估计方法 是一个值得关注的方向。
优点
能够提供较高的频率分辨率和较低的估计误差。
原理
格莱姆-梅尔谱估计利用了信号的模型参数，通过 构造一个模型函数来描述信号的频率响应特性， 并求解该函数的极值问题得到信号的功率谱。
缺点
需要预先设定模型函数的形式和参数，且计算复 杂度较高。
03 经典功率谱估计的优缺点
优点
01
02
03
算法成熟
经典功率谱估计方法经过 多年的研究和发展，已经 相当成熟，具有较高的稳 定性和可靠性。
9现代信号处理-功率谱估计(2)
rxx (1)
rxx(2)
rxx (0) rxx (1)
rxx
(N
1)
rxx
(N
2)
0
rxx
(
N
1)
rxx(N )
rxx (1)
可以看出AR模型得到的结果与按最大熵外推rxx(N+1)得到的结果 一致，这就证明了当x(n)为高斯分布时的最大熵谱估计与AR模型 法是等价的。
上式（4.6.8)是rxx(N+1)的一次函数，由此可解得rxx(N+1)。再 用类似的方法求得rxx(N+2)， rxx(N+3)，┄，然后确定功率谱估计。
最大熵谱估计用下式计算信号功率谱：
Pxx (e j )
2 w
N
1 ake jk
k 1
（4.6.9）
二、最大似然谱估计
１、最小方差谱估计
最大似然谱估计是用一个FIR滤波器实现，该滤波器对所
关心频率的正弦信号，可以无失真地通过，而对于其它频率的
信号，让其频响尽可能地小，亦即将它们尽可能地滤除。此时,
rxx (1)
rxx(2)
rxx (0) rxx (1)
rxx ( N
1)
rxx ( N
2) 0
（4.6.7）
rxx
(
N
1)
rxx(N )
rxx (1)
上 式 是 rxx(N+1) 的 一 次 函 数 ， 可 以 解 出 rxx(N+1) 。 继 续 再 将
rxx(N+1) 代 入 Rxx(N+2) 和 det(Rxx(N+2)) 中 ， 求 det(Rxx(N+2)) 对
现代信号处理大作业
现代信号处理大型作业一．试用奇阶互补法设计两带滤波器组(高、低通互补)，进而实现四带滤波器组；并画出其频响。
滤波器设计参数为：F p ＝1.7KHz , F r ＝2.3KHz , F s ＝8KHz , A rmin ≥70dB 。
（一）、分析与通常的滤波器相比，互补滤波器具有优良的结构特性和结构特性，具有较低的噪声能量和系数敏感性，其定义如下：一组滤波器H 12(),(),.......()Z H Z H Z n 如果满足下式：He Kjw k n(),==∑110<w<2π 则称这组滤波器为幅度互补滤波器;如果满足下式：He kjw k n()=∑=121, 0<w<2π则称这组滤波器为功率互补滤波器，同时互补滤波器还应该满足：Hz A z kk n()()=∑=1其中A(z)为全通函数，适当的选择全通函数，可以使两带函数具有所需要的低通和高通特性。
（二）、设计步骤(1) 对Fp 、Fr 进行预畸);();(''FsFrtg FsFptg r p ∏=Ω∏=Ω(2) 计算'''*r p c ΩΩ=Ω，判断'c Ω是否等于1，即该互补滤波器是否为互补镜像滤波器(3) 计算相关系数⎪⎩⎪⎨⎧-==+++=+-=-=ΩΩ=--=偶数)N 为(;21奇数)N 为 (;;lg /)16/1lg(;150152;1121;1;;])110)(110[(1213090500''02'''211－min1.0min1.0i i u q k N q q q q q k k q k k k k rp Ar Ap;)2cos()1(21))12(sin()1(21)1(21'2∑∑∞=∞=+-++-=Ωm mm m m m m i u Nm q u Nm q q ππ;42⎥⎦⎤⎢⎣⎡=N N;221N N N -⎥⎦⎤⎢⎣⎡=;)/1)(1(2'2'k k v i i i Ω-Ω-=12'1212,1;12N i v i i i =Ω+=--α 22'22,1;12N i v iii =Ω+=β (4) 互补镜像滤波器的数字实现;22i ii A αα+-=;22iii B ββ+-=1221,1;1)(N i ZA Z A Z H i i i =++=∏--22212,1;1)(N i ZB Z B Z Z H i i i =++=∏--- )];()([21)(21Z H Z H Z H L +=（三）、程序与结果 1． 二带滤波器组 （1） 源程序： clear; clf;Fp=1700;Fr=2300;Fs=8000; Wp=tan(pi*Fp/Fs); Wr=tan(pi*Fr/Fs); Wc=sqrt(Wp*Wr); k=Wp/Wr;k1=sqrt(sqrt(1-k^2)); q0=0.5*(1-k1)/(1+k1);q=q0+2*q0^5+15*q0^9+150*q0^13; N=11;N2=fix(N/4); M=fix(N/2); N1=M-N2; for jj=1:M a=0;for m=0:5a=a+(-1)^m*q^(m*(m+1))*sin((2*m+1)*pi*jj/N);%N is odd, u=j end ab=0;for m=1:5b=b+(-1)^m*q^(m^2)*cos(2*m*pi*jj/N); end bW(jj)=2*q^0.25*a/(1+2*b);V(jj)=sqrt((1-k*W(jj)^2)*(1-W(jj)^2/k)); endfor i=1:N1alpha(i)=2*V(2*i-1)/(1+W(2*i-1)^2); endfor i=1:N2beta(i)=2*V(2*i)/(1+W(2*i)^2); endfor i=1:N1a(i)=(1-alpha(i)*Wc+Wc^2)/(1+alpha(i)*Wc+Wc^2); endfor i=1:N2b(i)=(1-beta(i)*Wc+Wc^2)/(1+beta(i)*Wc+Wc^2); endw=0:0.0001:0.5;LP=zeros(size(w));HP=zeros(size(w));for n=1:length(w)z=exp(j*w(n)*2*pi);H1=1;for i=1:N1H1=H1*(a(i)+z^(-2))/(1+a(i)*z^(-2)) ;endH2=1/z;for i=1:N2H2=H2*(b(i)+z^(-2))/(1+b(i)*z^(-2));endLP(n)=abs((H1+H2)/2);HP(n)=abs((H1-H2)/2);endplot(w,LP,'b',w,HP,'r');hold on;xlabel('digital frequency');ylabel('amptitude');（2）运行结果：见图1图1 二带数字滤波器组2．四带滤波器组（1）源程序：clf;Fp=1700;Fr=2300;Fs=8000;Wp=tan(pi*Fp/Fs);Wr=tan(pi*Fr/Fs);Wc=sqrt(Wp*Wr);k=Wp/Wr;k1=sqrt(sqrt(1-k^2));q0=0.5*(1-k1)/(1+k1);q=q0+2*q0^5+15*q0^9+150*q0^13;N=11;N2=fix(N/4);M=fix(N/2);N1=M-N2;for jj=1:Ma=0;for m=0:5a=a+(-1)^m*q^(m*(m+1))*sin((2*m+1)*pi*jj/N); % N is odd, u=jendb=0;for m=1:5b=b+(-1)^m*q^(m^2)*cos(2*m*pi*jj/N);endW(jj)=2*q^0.25*a/(1+2*b);V(jj)=sqrt((1-k*W(jj)^2)*(1-W(jj)^2/k));Endfor i=1:N1alpha(i)=2*V(2*i-1)/(1+W(2*i-1)^2);endfor i=1:N2beta(i)=2*V(2*i)/(1+W(2*i)^2);endfor i=1:N1a(i)=(1-alpha(i)*Wc+Wc^2)/(1+alpha(i)*Wc+Wc^2);endfor i=1:N2b(i)=(1-beta(i)*Wc+Wc^2)/(1+beta(i)*Wc+Wc^2);endw=0:0.0001:0.5;LLP=zeros(size(w));LHP=zeros(size(w));HLP=zeros(size(w));HHP=zeros(size(w));for n=1:length(w)z=exp(j*w(n)*2*pi);H1=1;for i=1:N1H1=H1*(a(i)+z^(-2))/(1+a(i)*z^(-2)) ;endH21=1;for i=1:N1H21=H21*(a(i)+z^(-4))/(1+a(i)*z^(-4)) ;H2=1/z;for i=1:N2H2=H2*(b(i)+z^(-2))/(1+b(i)*z^(-2));endH22=1/(z^2);for i=1:N2H22=H22*(b(i)+z^(-4))/(1+b(i)*z^(-4));endLP=((H1+H2)/2);HP=((H1-H2)/2);LLP(n)=abs((H21+H22)/2*LP);LHP(n)=abs((H21-H22)/2*LP);HHP(n)=abs((H21+H22)/2*HP);HLP(n)=abs((H21-H22)/2*HP);endplot(w,LLP,'b',w,LHP,'r',w,HLP,'k',w,HHP,'m')hold onxlabel('digital frequency');ylabel('amptitude');（2）运行结果：见图2图2 四带数字滤波器组二、根据《现代数字信号处理》第四章提供的数据，试用如下方法估计其功率谱，并画出不同参数情况下的功率谱曲线：1）Levison算法2）Burg算法3） ARMA 模型法 4） MUSIC 算法 1 Levinson 算法Levinson 算法用于求解Yule-Walker 方程，是一种按阶次进行递推的算法，即首先以AR （0）和AR （1）模型参数作为初始条件，计算AR （2）模型参数；然后根据这些参数计算AR （3）参数，等等，一直到计算出AR （p ）模型参数为止，需要的运算量数量级为2p ，其中p 为AR 模型的阶数。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
清华大学《数字信号处理》期末作业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 xk SN R y k x k =-∑∑其中，()x k 为原nosidopp 信号，()y k 为滤波后信号。
SNR 越大表示滤除部分能力越小，可以反映滤波后信号对原信号的跟踪能力，对前128点主要考察SNR ，越大滤波器性能越好。
对128点以后的点，信号能量集中在低频部分，滤波后效果不仅仅需考察SNR ，同时考察滤波后信号的平滑程度。
定义平滑指标[1]22((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 Mx x x Mx L x L x N-⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥--⎣⎦X得到的矩阵为方阵可最大程度利用其特征值，对nosidopp信号，N=1024，故选择L=513，M=512构造矩阵。
利用svd分解得到特征值，如图1.21。
同时得到矩阵U，V。
图1.22 矩阵特征值可以看出特征值存在明显的转折点，将转折点（第30点）以后的特征值值为0。
恢复信号矩阵：1112111111ˆˆˆ(0)(1)(1)ˆˆˆ(1)(2)()ˆˆˆˆ(1)()(1)T xx xM x x xM xL xL xN -⎡⎤⎢⎥⎢⎥==⎢⎥⎢⎥--⎣⎦X U SV此时得到的矩阵1ˆX 将不再对称，按上式中对应位置得到滤波后信号ˆ()y x n =为：1ˆˆ()(),0,1, (1)kxi 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 Mx x x Mx 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)。