现代信号分析课程大作业 - 图文

课程名称： 现 代 信 号 分 析 考试形式：□专题研究报告 □论文 大作业 □综合考试 学生姓名： 学号：

序 号 1 2 3 4 5 6 7 8 9 10 分 项 类 别 题目1 题目2 题目3 题目4 题目5 题目6 题目7 题目8 题目9 题目10 得 分 评阅人： 时 间： 年 月

总 分

目录

一、FFT算法的特点。 ......................................... - 3 - 二、平稳随机过程 ............................................. - 8 - 三、功率谱估计 .............................................. - 10 - 四、自适应滤波器（LMS和RLS） ............................... - 16 - 五、维纳滤波器 .............................................. - 23 - 六、FIR维纳滤波器 .......................................... - 25 - 七、AR/MA滤波器 ............................................ - 28 - 八、卡尔曼滤波器 ............................................ - 31 - 九、小波分析方法及应用 ...................................... - 37 - 十、其他现代信号分析方法 .................................... - 43 -

现代信号分析

一、FFT算法的特点。

题目：请举例说明FFT算法的特点

计算离散傅里叶变换的一种快速算法，简称FFT。快速傅里叶变换是1965年由J.W.库利和T.W.图基提出的。采用这种算法能使计算机计算离散傅里叶变换所需要的乘法次数大为减少，特别是被变换的抽样点数N越多，FFT算法计算量的节省就越显著。根据对序列分解与选取方法的不同而产生了FFT的多种算法，基本算法是基2DIT和基2DIF。FFT在离散傅里叶反变换、线性卷积和线性相关等方面也有重要应用。

快速傅氏变换（FFT），是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。有限列长为N的序列x（n）的DFT变换公式在MATLAB中的表达式为

nkX(k)??x(n)WNn?1N?1

k?1,2,...,N?1

其逆变换为

1N?nkx(n)??X(k)WN n?1,2,...,N

Nk?1FFT主要就是利用了WN以下两个特性使长序列的DFT分解为更小点数的DFT来实现的。 （1）

利用WN的对称性使DFT运算中有些项合并

k(N?n)?knkn*WN?WN?(WN)

nknk（2）

利用WN的周期性和对称性使长序列的DFT分解为更小点数的DFT

nkk(n?N)(k?N)nWN?WN?WN

nkFFT算法在应用时能够大大减少计算量，但是也存在着和DFT同样的一些问题，主要有：频谱混

叠，泄露问题以及栅栏效应的问题，下面通过运用MATLAB的FFT相关函数对信号进行变换来进行实际分析。

（1） 频谱混叠

奈奎斯特定理已被众所周知了，所以几乎所有人的都知道为了不让频谱混叠，理论上采样频谱大于等于信号的最高频率。采样周期的倒数是频谱分辨率，最高频率的倒数是采样周期。设定采样点数为N，采样频率fs，最高频率fh，故频谱分辨率f=fs/N,而fs>=2fh，所以可以看出最高频率与频谱分辨率是相互矛盾的，提高频谱分辨率f的同时，在N确定的情况下必定会导致最高频率fh的减小；同样的，提高最高频率fh的同时必会引起f的增大，即分辨率变大,下面举例说明。

设输入信号为s?3cos240?t?5sin300?t，即它的主要频率为120Hz和150Hz。分别取采样频率为400Hz，300Hz，200Hz，100Hz，程序如下： L=1024; %采样点数 NFFT = 2^nextpow2(L); % Next power of 2 from length of y Fs1=400; % 采样频率 - 3 -

t1=(0:L-1)/Fs1; s1=3*cos(2*pi*120*t1)+5*sin(2*pi*150*t1); Y1 = fft(s1,NFFT)/L; f1 = Fs1/2*linspace(0,1,NFFT/2+1); Fs2=300; t2=(0:L-1)/Fs2; s2=3*cos(2*pi*120*t2)+5*sin(2*pi*150*t2); Y2= fft(s2,NFFT)/L; f2=Fs2/2*linspace(0,1,NFFT/2+1); Fs3=200; t3=(0:L-1)/Fs3; s3=3*cos(2*pi*120*t3)+5*sin(2*pi*150*t3); Y3 = fft(s3,NFFT)/L; f3=Fs3/2*linspace(0,1,NFFT/2+1); Fs4=100; t4=(0:L-1)/Fs4; s4=3*cos(2*pi*120*t4)+5*sin(2*pi*150*t4); Y3 = fft(s4,NFFT)/L; f4=Fs4/2*linspace(0,1,NFFT/2+1); subplot(2,2,1); plot(f1,2*abs(Y1(1:NFFT/2+1))); title('采样频率为400Hz'); grid on; xlabel('频率(Hz)'); subplot(2,2,2); plot(f2,2*abs(Y2(1:NFFT/2+1))); title('采样频率为300Hz'); grid on; xlabel('频率(Hz)'); subplot(2,2,3); plot(f3,2*abs(Y3(1:NFFT/2+1))); title('采样频率为200Hz'); grid on; xlabel('频率(Hz)'); subplot(2,2,4); plot(f4,2*abs(Y4(1:NFFT/2+1))); title('采样频率为100Hz'); grid on; xlabel('频率(Hz)'); grid on; 结果如图1-1所示，可见当采样频率为300Hz时，没有发生频率混叠，但当频率小于两倍的信号频率时，经过采样就会发生混叠现象，输入信号中的某些频率就无法辨识出。

- 4 -

现代信号分析

图1-1 FFT的混淆现象

（2） 泄漏现象

实际中的信号序列往往很长，甚至是无限长序列。为了方便，我们往往用截短的序列来近似它们。这样可以使用较短的DFT来对信号进行频谱分析。这种截短等价于给原信号序列乘以一个窗函数。而窗函数的频谱不是有限带宽的，从而它和原信号的频谱进行卷积以后会扩展原信号的频谱。值得一提的是，泄漏是不能和混叠完全分离开的，因为泄露导致频谱的扩展，从而造成混叠。为了减小泄漏的影响，可以选择适当的窗函数使频谱的扩散减到最小。下面举例说明。

设输入信号为s?3cos240?t?5sin300?t，即它的主要频率为120Hz和150Hz。，用MATLAB自带的信号处理工具箱（通过命令sptool调出）分别对信号加矩形窗（boxcar），海明窗（Hamming），汉宁窗（Hanning）和布莱克曼窗（Blackman），分析不同的窗函数对减小泄露现象的影响。

具体操作如下：先生成信号s，然后在sptool工具箱的file菜单中import信号s，Fs为400Hz，如图1-2所示。

图1-2 sptool工具箱界面

然后点create按钮，在patameters栏选择welch方法，Nfft点数1024，Nwind为256，分别加布

- 5 -

莱克曼窗,矩形窗，海明窗，和三角窗，得如下PSD图。

(a)blackman (b)boxcar

(c)hamming (d)triang

图1-3 FFT的泄露现象

由以上四图可见，布莱克曼窗主瓣宽，旁瓣小，频率识别精度最低，但幅值识别精度最高。矩形窗主瓣窄，旁瓣大，频率识别精度最高，幅值识别精度最低；海明窗是余弦窗的一种，又称改进的升余弦窗，海明窗与汉宁窗都是余弦窗，只是加权系数不同。海明窗加权的系数能使旁瓣达到更小，海明窗的频谱是由3个矩形时窗的频谱合成，但其旁瓣衰减速度为20dB／（10oct）。三角窗与矩形窗比较,主瓣宽约等于矩形窗的两倍,但旁瓣小,而且无负旁瓣。

所以对于窗函数的选择，应考虑被分析信号的性质与处理要求。如果仅要求精确读出主瓣频率，而不考虑幅值精度，则可选用主瓣宽度比较窄而便于分辨的矩形窗，如果分析窄带信号，且有较强的干扰噪声，则应选用旁瓣幅度小的窗函数，如汉宁窗。

（3） 栅栏效应

栅栏效应是因为用计算频谱只限制为基频的整数倍而不可能将频谱视为一连续函数而产生的。设信号的最高频率为

fh，要求的频率分辨率为F（即在频率轴上所能得到的最小频率间隔），则采样点数N需要满足以下条件：

N?2fhF

若不满足，就会产生栅栏效应。减小栅栏效应通常有两个方法：

一个方法是在原纪录末端添加一些零值点来变动时间周期内的点数，并保持记录不变。从而在保持原有频谱连续形式不变的情况下，变更了谱线的位置。这样，原来看不到的频谱分量就能移到可见的位置上了。但是，补零并没有增加序列的有效长度，所以并不能提高分辨率。但补零同时可以使数据长度N为2的整数幂，以便于应用FFT。

另一个方法是增加采样点数，从而减小频率间隔，使原来被挡住的一些频谱的峰点或谷点显露出来。注意，这时候每根谱线多对应的频率和原来的已经不相同了。 下面通过实际的函数来分析两种方法的效果。

- 6 -

现代信号分析

程序如下：

L=90; % 采样点数 NFFT = 2^nextpow2(L); % Next power of 2 from length of y Fs=400; % 2采样频率 t=(0:L-1)/Fs; s=cos(200*pi*t)+2*cos(204*pi*t); Y1 = fft(s,NFFT)/L; f1 = Fs/2*linspace(0,1,NFFT/2+1); plot(f1,2*abs(Y1(1:NFFT/2+1))); xlabel('频率'); % ----原序列末端补零 ------------------------------------------- L3=1024; NFFT= 2^nextpow2(L3); s3=[s zeros(1,L3-L)]; Y3= fft(s3,NFFT)/L3; f3=Fs/2*linspace(0,1,NFFT/2+1); figure; plot(f3,2*abs(Y3(1:NFFT/2+1))); xlabel('频率'); ('频率'); % ---- 增加采样点数--------------------------------------------- L2=1024; NFFT = 2^nextpow2(L2); t=(0:L2-1)/Fs; s2=cos(200*pi*t)+2*cos(204*pi*t); Y2 = fft(s2,NFFT)/L2; f2 = Fs/2*linspace(0,1,NFFT/2+1); figure; plot(f2,2*abs(Y2(1:NFFT/2+1))); xlabel('频率'); 设输入信号为s?cos200?t?2cos204?t，采样频率为400Hz。可见，信号含有2个频率分量，分别为100Hz和102Hz，则要求的频率分辨率为2Hz，可算得N为202，即至少采样202个点，才不会发生栅栏现象。

首先取采样点数为90，对其进行FFT，得到的图形如图1-4所示。

图1-4 FFT的栅栏效应

- 7 -

可见，虽然采样频率满足香农采样定理，但由于采样点数较少，频率为20Hz和20.5Hz的信号无法分辨。下面采用补零的方法，将N补满1024，得到的FFT图形如图1-5所示。

图1-5 信号末端补零改善栅栏效应

由上图可以看出，补零对分辨率没有影响，但是对频谱起到了平滑的作用。同时，由于增加了大量零值，使得信号的平均幅频响应大大减小。

然后采用增加采样点的方法，将采样点数增加到1024，FFT后的图形如图1-6所示。

图1-6 增加采样点数改善栅栏效应

由上图可以看出，增加采样点数后就能分辨出100Hz和102Hz的信号，有效地改善了栅栏效应。从时域看，这种方法相当于对信号进行整周期采样，实际中常用此方法来提高周期信号的频谱分析精度。

二、平稳随机过程

题目：什么是宽平稳随机过程？什么是严平稳随机过程？它们之间有什么联系？（5分）

平稳随机过程是指它的统计特性不随时间的推移而变化，在谱估计中可以认为平稳随机信号是由一个白噪声激励一个线性时不变系统产生的。

- 8 -

现代信号分析

2.1宽平稳随机过程

设有一个二阶随机过程{?(t),t?T}，它的均值为常数，自相关函数仅是时间差的函数，则称它为宽平稳随机过程，用符号化语言表示，给定一个二阶矩过程{?(t),t?T}，如果对任意t,t???T，有：

（1）E[?(t)]?a（常数）

（2）自相关函数R(t,t??)?E[?(t)?(t??)]?R(?) 则称{?(t),t?T}为宽平稳随机过程或广义平稳随机过程。

2.2 严平稳随机过程

在数学中，平稳随机过程（Stationary random process）或者严平稳随机过程（Strictly-sense stationary random process），又称狭义平稳随机过程，是指随机过程在某一固定时间和位置的概率分布与所有时间和位置的概率分布相同，即随机过程的统计特性不随时间的推移而变化。这样，数学期望和方差这些参数也不随时间和位置变化。

用符号化语言表示出来，设随机过程?，对于任意的n（n=1,2,···）和任意选定的?以及?为任意值，当?时，n维随机变量?和??(t1??),?(t2??),???,?(tn??)?具有相同的概率分布，即：

fn(x1,x2,?,xn;t1,t2,?,tn)?fn(x1,x2,?,xn;t1??,t2??,?,tn??)

则称随机过程??(t),t?T?具有平稳性，称此过程为严平稳随机过程，简称随机过程。

2.3 严平稳随机过程与宽平稳随机过程区别联系

（1）一个宽平稳过程不一定是严平稳过程，一个严平稳过程也不一定是宽平稳过程。 例1：X(n)?sinwn,n?0,1,2???，其中w服从U(0,2π)，随机过程?X(n),n?0,1,2????是宽平稳过程，但不是严平稳过程；

例2：服从柯西分布的随机变量序列是严平稳随机过程，但不是宽平稳随机过程。

（2）宽平稳过程定只涉及与一维、二维分布有关的数字特征，所以一个严平稳过程只要二阶矩存在（或E[?(t)]有界），则必定是宽平稳随机过程。但反过来一般不成立。

2- 9 -

注：二阶矩过程定义，如果随机过程??(t),t?T?对第一个t?T，二阶矩E[?(t)??(t)]都存在，那么称其为二阶矩过程。

（3）正态过程是一个重要特例，一个宽平稳的正态过程必定是严平稳的。这是因为：正态过程的概率密度是由均值函数和自相关函数完全确定的，因而如果均值函数和自相关函数不随时间的推移而变化，则概率密度函数也不随时间的推移发生变化。

三、功率谱估计

题目：简述经典功率谱估计与现代功率谱估计的差别，并计算下题：（14分）

设随机过程是单位方差白噪声激励如下的系统而产生的：

s(n)?2.7377s(n?1)?3.7476s(n?2)?2.6293s(n?3)?0.9224s(n?4)??(n)

取序列的长度分别为128和1024数据段，并用经典和现代谱估计方法进行谱估计，分析不同参数对最终结果的影响。 解答：

3.1 经典功率谱估计和现代功率谱估计的差别

信号的频谱分析是研究信号特性的重要手段之一，通常是求其功率谱来进行频谱分析。功率谱反映了随机信号各频率成份功率能量的分布情况，可以提示信号中隐含的周期性及靠得很近的谱峰等有用信息，在许多领域都发挥了重要作用。然而，实际应用中的平稳随机信号通常是有限长的，只能根据有限长信号估计原信号的真实功率谱。因此，功率谱估计就是基于有限的数据寻找信号、随机过程或系统的频率成分。一般功率谱估计分为经典谱估计和现代谱估计。其各自的特点及相互差别如下。

根据所用理论的不同，通常将基于相关函数的傅里叶变换的估计方法称为经典功率谱估计，而将参数模型估计方法和基于相关矩阵特征值分解的信号频率估计方法，称为现代功率谱估计方法。经典谱估计为线性估计方法，是建立在传统的傅里叶变换的基础之上，它首先由给定的数据估计自相关序列rx(k)，然后对估计出的自相关序列rx(k)进行傅里叶变换得到功率谱估计，没有将信号的可用信息结合到估计过程中，有计算效率高、估计值正比于正弦波信号的功率等优势。但性能依赖于数据序列的长度，频率分辨率低，且需在方差和分辨率之间做出权衡，不适用于短时数据的情况，典型代表有Blackman和Tukey提出的自相关谱估计法（简称BT法），周期图法以及周期图法的改进方法——Barlett法和Welth法；而现代功率谱估计是非线性估计方法，估计性能依赖于参数模型，最终可获得方差小、分辨率高的谱估计，但所用模型必须适合于所分析信号，否则谱估计将是错误或者不准确的。有AR模型，MA模型，ARMA模型等；非参数模型谱估计有最小方差法和MUSIC法等。

- 10 -

?

现代信号分析

3.2功率谱计算

由上式的输入输出关系可以看出这是一个AR模型，由已知数据得到：

H(z)?11?2.2377z?1?3.7476z?2?2.6293z?3?0.9224z?4

本题中的随机过程是由单位方差白噪声激励一个差分系统方程而产生。首先利用Matlab产生该信号源：

n=1:1:1024; wn=randn(1,length(n)); %产生随机单位白噪声序列 sn=filter([1],[1,-2.2377,3.7476,-2.6293,0.9224],wn); subplot; plot(sn,'r','LineWidth',1.4); title('s(n)','fontsize',12); 产生的信号源如下图3-1所示。对信号源信号求FFT变换，得到信号源的频谱图，作为估计算法的参考标准，如图3-2所示。

图3-1 Matlab产生的信号源 图3-2 信号源的频谱图

1）经典谱估计

下面首先采用直接法（周期图法）进行功率谱估计，程序如下： % nfft=128 window=boxcar(length(sn)); nfft=128; [Pxx,f]=periodogram(sn,window,nfft); subplot(211);plot(f,10*log10(Pxx)); xlabel('频率');ylabel('相对功率谱密度dB'); title('周期图法 点数128'); % nfft=1024 window=boxcar(length(sn)); % 矩形窗 nfft=1024; [Pxx,f]=periodogram(sn,window,nfft); subplot(212);;plot(f,10*log10(Pxx)); xlabel('频率'); ylabel('相对功率谱密度dB'); title('周期图法 点数1024'); 改变采样点数，获得的周期图法谱估计如下图所示： - 11 -

图3-3 不同序列长度周期法谱估计

运用周期图法通过对采样数据长度的调整，即在采样数据为128、1024时比较周期图法的效果，可知采集数据序列长度越大，谱分辨率越高，但谱曲线起伏也越大。因此需要改进，这里的改进指的是数据方差的改进。当前主要有两种改进方法：一种是周期图的平滑，即采用间接法估计功率谱，即BT法；另一种是平均法，它将长度为N的数据X(n)分成L段，分别求出每段的功率谱，然后加以平均，如Barlett法和Welch法。这里我们就采用以上3种方法对直接法进行改进。Matlab程序如下：

%BT法 nfft=1024; %以点数1024为例 cxn=xcorr(sn,'unbiased'); %求序列自相关函数 P1=fft(cxn,nfft); subplot(311); plot(f,10*log10(P1(1:257)),'r','LineWidth',1.4); xlabel('频率','fontsize',12);ylabel('功率谱/dB ','fontsize',12); title('BT法','fontsize',12) % Bartlett法 window=boxcar(length(sn)); %矩形窗 noverlap=0; %数据无重叠 p=0.9; %置信概率 p2=psd(sn,nfft,1,window,noverlap,p); %计算PSD subplot(312);plot(f,10*log10(p2(1:257)),'r','LineWidth',1.4); xlabel('频率','fontsize',12);ylabel('功率谱/dB','fontsize',12); title('Bartlett法','fontsize',12); %Welch法 window=boxcar(length(sn)); %矩形窗 noverlap=0; %数据无重叠 range='onesided'; %频率间隔为[0,1/2]，只计算一半频率 [p3,w]=pwelch(sn,window,noverlap,nfft,range); subplot(313);plot(f,10*log10(p3),'b','LineWidth',1.4); xlabel('频率','fontsize',12);ylabel('功率谱/dB','fontsize',12); title('Welch法','fontsize',12); 各个经典功率谱估计方法的结果如图3-4所示。 - 12 -

现代信号分析

图3-4 经典谱估计方法的估计结果

为了了解使用不同种窗函数对功率谱的影响，以Welth法为例，进行了应用各种常见窗函数的仿真实验。程序如下:

%矩形窗 window=boxcar(length(sn)); %矩形窗 noverlap=0; %数据无重叠 range='onesided'; %频率间隔为[0,1/2]，只计算一半频率 [p3,w]=pwelch(sn,window,noverlap,nfft,range); subplot(221);plot(f,10*log10(p3),'b','LineWidth',1.4); xlabel('频率','fontsize',12);ylabel('功率谱/dB','fontsize',12); title('矩形窗','fontsize',12); %汉明窗 window=hamming(length(sn)); %矩形窗 noverlap=0; %数据无重叠 range='onesided'; %频率间隔为[0,1/2]，只计算一半频率 [p3,w]=pwelch(sn,window,noverlap,nfft,range); subplot(222);plot(f,10*log10(p3),'b','LineWidth',1.4); xlabel('频率','fontsize',12);ylabel('功率谱/dB','fontsize',12); title('汉明窗','fontsize',12); %布莱克曼窗 window=blackman(length(sn)); %矩形窗 noverlap=0; %数据无重叠 range='onesided'; %频率间隔为[0,1/2]，只计算一半频率 [p3,w]=pwelch(sn,window,noverlap,nfft,range); subplot(223);plot(f,10*log10(p3),'b','LineWidth',1.4); xlabel('频率','fontsize',12);ylabel('功率谱/dB','fontsize',12); title('布莱克曼窗','fontsize',12); %三角窗 window=triang(length(sn)); %矩形窗 noverlap=0; %数据无重叠 range='onesided'; %频率间隔为[0,1/2]，只计算一半频率 [p3,w]=pwelch(sn,window,noverlap,nfft,range); subplot(224);plot(f,10*log10(p3),'b','LineWidth',1.4); - 13 -

xlabel('频率','fontsize',12);ylabel('功率谱/dB','fontsize',12); title('三角窗','fontsize',12); 估计结果如图3-5所示：

图3-5 不同窗函数对功率谱估计结果的影响

综合以上的波形，可以对经典功率谱估计作如下总结：经典功率谱估计，不论是周期图法还是BT法都可以用FFT快速计算，且计算量较小，物理概念明确，因而仍是目前较常用的谱估计方法。功率谱的分辨率正比于数据长度，当数据长度较大时，谱的曲线起伏加剧。BT法是对周期图法的改进，但其谱的平滑是以牺牲分辨率为代价的。Welth方法是对BT法德改进，它能使用更多的窗函数，并能更好的控制功率谱密度的估计的方差特性。在实际应用中，在方差、偏差和分辨率之间存在着矛盾，只能根据需要作折中处理。

2）现代谱估计

下面，根据AR模型法和burg算法对信号源进行现代功率谱估计。Matlab程序如下：

%------------现代功率谱分析---------% n=1:1:1024; wn=randn(1,length(n)); %产生随机单位白噪声序列 sn=filter([1],[1,-2.7377,3.7476,-2.6293,0.9224],wn); [Pxx,f]=periodogram(sn,window,nfft); nfft=1024; %AR模型法 [p4,w]=pyulear(sn,100,nfft); figure(3);subplot(2,1,1); plot(f,10*log10(abs(p4)),'b','LineWidth',1.4); xlabel('频率 ','fontsize',12); ylabel('功率谱/dB','fontsize',12); title('AR模型法','fontsize',12); %burg算法 [p5,w]=pburg(sn,100,nfft); subplot(2,1,2); plot(f,10*log10(abs(p5)),'b','LineWidth',1.4); xlabel('频率','fontsize',12); ylabel('功率谱/dB','fontsize',12); title('burg算法','fontsize',12); - 14 -

现代信号分析

估计结果如图3-6所示：

图3-6 现代功率谱估计方法的估计结果

为了分析现代功率谱中模型阶数对谱估计的影响，现利用AR模型法，改变模型阶数。分别用10、20、100作为模型阶数在采样点数分别为128和1024下进行分析，结果如下： %---------=-------现代功率谱分析---=---------% %AR模型法 n=1:1:1024; wn=randn(1,length(n)); %产生随机单位白噪声序列 sn=filter([1],[ 1,-2.7377,3.7476,-2.6293,0.9224],wn); [P1,w4]=pmcov(sn,10,128); [P2,w5]=pmcov(sn,50,128); [P3,w6]=pmcov(sn,100,128); [P4,w1]=pmcov(sn,10,1024); [P5,w2]=pmcov(sn,50,1024); [P6,w3]=pmcov(sn,100,1024); subplot(3,2,1); plot(w4/pi,10*log10(P1)); xlabel('频率 ','fontsize',12); ylabel('功率谱/dB','fontsize',12); title('128 10阶AR模型功率谱图');grid on; subplot(3,2,3); plot(w5/pi,10*log10(P2)); xlabel('频率 ','fontsize',12); ylabel('功率谱/dB','fontsize',12); title('128 50阶AR模型功率谱图');grid on; subplot(3,2,5); plot(w6/pi,10*log10(P3)); xlabel('频率 ','fontsize',12); ylabel('功率谱/dB','fontsize',12); title('128 100阶AR模型功率谱图');grid on; subplot(3,2,2); plot(w1/pi,10*log10(P4)); xlabel('频率 ','fontsize',12); ylabel('功率谱/dB','fontsize',12); title('1024 10阶AR模型功率谱图');grid on; subplot(3,2,4); plot(w2/pi,10*log10(P5)); xlabel('频率 ','fontsize',12); ylabel('功率谱/dB','fontsize',12); title('1024 50阶AR模型功率谱图');grid on; subplot(3,2,6); plot(w3/pi,10*log10(P6)); xlabel('频率 ','fontsize',12); ylabel('功率谱/dB','fontsize',12); title('1024 100阶AR模型功率谱图');grid on; - 15 -

图3-7 不同阶数下AR模型法进行功率谱估计

可以看到，在相同采样数据长度下阶数过小会导致频率成分的辨识出现偏差，阶数过大会使随机起伏变化密集和加剧，方差增大。在50阶的情况下效果较为理想，但阶数越高的情况下计算量同时也在增加，应该适当的进行选择。在相同AR模型阶数的条件下，采样数据越多，频率的分辨率也越高，但会出现波动。

由上述分析可以看出，AR模型法、burg算法估计出的功率谱相对经典的功率谱估计方法得出的结果方差较小，频率的准确性也相对较高。现代谱估计的分辨率可以不受样本长度的限制。这是因为对于给定的N 点有限长序列x(n), 虽然其估计出的自相关函数也是有限长的, 但是现代谱估计的一些隐含着数据和自相关函数的外推, 使其可能的长度超过给定的长度,不像经典谱估计那样受窗函数的影响。与经典谱估计结果比较，现代谱估计在分辨率和方差上都优于经典谱，而且现代谱线要平滑得多，且在采样数据短的情况下也非常适用。

四、自适应滤波器（LMS和RLS）

题目：请简述一下自适应滤波器（LMS和RLS），并任选下面的一题，（2）、（3）数据可从网上下载（MIT-BIH

数据库：

http://www.physionet.org/physiobank/database/mitdb），分析步长及指数加权因子对滤波结果的影响。同时给出迭代次数与滤波器系数，迭代次数与均方误差之间的关系曲线。（13分）

（1）利用自适应滤波器从宽带信号中提取单频信号。

s(t)是宽频信号， 设x(t)?s(t)?cos(2?ft)x(t)?s(t)?Acos(2?f1t??)?Bcos(2?f2t)，

A,B, f1, f2任选，要求提取两个单频信号；同时对不同的信噪比的情况比较一下。

（2）消除心电图的电源干扰 （3）胎儿心电监护

- 16 -

现代信号分析

滤波技术是指从复杂的含有噪声的信号中提取出有用的所需信号，而滤波器是一个选频系统，让有用的所需信号以尽量小的衰减通过，其实我们可以把大部分的信号处理系统看作滤波器系统。滤波器可以分为经典滤波器和现代滤波器，也可以分为模拟滤波器和数字滤波器，而自适应滤波器是属于现代滤波器，它是相对固定滤波器而言。本题首先介绍自适应滤波器的概念和工作原理，然后对LMS和RLS两种算法进行介绍，随后以胎儿心电信号（FECG）的提取（题3）为例来介绍自适应滤波器的应用，并分析步长及指数加权因子对滤波结果的影响。 1. 自适应滤波器简介

自适应滤波器是指能够根据输入信号自动调整性能进行数字信号处理的数字滤波器。自适应滤波器在许多领域都有应用，如自适应天线系统、数字通信接收机、自适应噪声对消技术、系统建模等，本题的例子便是应用自适应滤波器进行自适应噪声对消技术。

自适应滤波器通常由两个部分组成，一部分是数字滤波器，用来完成希望的信号处理；另一部分是自适应算法，用来调节滤波器系数（或加权），如图3-1所示。

d(n) + ∑ - x(n) 参数可调数字滤波器 自适应滤波器算法 e(n) y(n) 图4-1 自适应滤波器的结构框图

如图，参数可调数字滤波器可以是FIR数字滤波器或IIR数字滤波器，也可以是格型数字滤波器。输入信号x(n)通过参数可调数字滤波器后产生输出信号(或响应)y(n)，将其与参考信号(或称期望响应)d(n)进行比较，形成误差信号e(n)。e(n)(有时还要利用x(n))通过某种自适应算法对滤波器参数进行调整，最终使e(n)的均方值最小。根据算法的不同，可以将自适应滤波器分为：递推最小二乘（RLS）滤波器、最小均方差（LMS）滤波器、格型滤波器、无限冲激响应（IIR）滤波器等。

（1）LMS滤波器

最小均方误差(LMS)算法是由Widrow和Hoff提出的，具有计算量小、易于实现等优点在实践中广泛应用。LMS算法的基本思想：调整滤波器自身的参数，使滤波器的输出信号与期望输出信号之间的均方误差最小，系统输出为有用信号的最佳估计。实质上LMS可以看成是一种随机梯度或者随机逼近算法，可以写成如下的基本迭代方程：

- 17 -

M?1?Ty(k)?W(k)X(k?i)?i?i?0?? ?e(k)?d(k)?y(k)?W(k?1)?W(k)?2?e(k)X(k?i)i?i??其中μ为步长因子，是控制稳定性和收敛速度的参量，X(k)表示第k时刻参考信号矢量，

X(k)?[n(k),n(k?1)???n(k?M?1)]，k为迭代次数，M为滤波器的阶数。d(k)表示第k时刻

的输入信号矢量，y(k)、e(k)分别表示第k时刻的输出信号与输出误差，W(k)表示k时刻权系数矢量，W(k)=[W(k，0)，W(k，1)…W(k，M-1)]。从上式可以看出其结构简单、计算量小且稳定性好等优点，但固定步长的LMS算法在收敛速度、跟踪速率及权失调噪声之间的要求是相互矛盾的，为了克服其缺点，人们提出了各种变步长的LMS改进算法，主要是采用减小均方误差或者以某种规则基于时变步长因子跟踪信号的时变，其中有正规LMS算法(NLMS)、梯度自适应步长算法、自动增益控制自适应算法、符号-误差LMS算法、符号-数据LMS算法、数据复用LMS算法等。

自适应滤波器收敛的条件是0???1?max。其中?max是输入信号的自相关矩阵R的最大特征值。

μ的选取必须在收敛速度和失调之间取得较好的折中，既要具有较快的收敛速度，又要使稳态误差最小。它控制了算法稳定性和自适应速度，如果μ很小，算法的自适应速度会很慢；如果μ很大，算法会变得不稳定。

（2）RLS滤波器

RLS（Recursive Least-Square）算法是通过矩阵求逆引理而导出的，它对输入信号的自相关矩阵

Rxx(n)的逆矩阵进行递推估计更新，其迭代公式如下：

?e(n)?d(n)?xT(n)w?(n?1),?P(n?1)x(n)?g(n)?,?T??x(n)P(n?1)x(n) ??P(n)???1P(n?1)???1g(n)xT(n)P(n?1),???(n)?w?(n?1)?g(n)e(n)?w上式中，n是递推次数，输入为x(n)，P(n)为信号自相关矩阵的逆，g(n)为增益矢量，误差为e(n)，期望输出为d(n)，权重因子为λ。

RLS算法能实现快速收敛，且其收敛速率不随输入向量x ( n)集平均相关矩阵R特征值的扩散度而改变。当工作于时变环境中时，这类算法具有极好的性能，但其实现都以增加计算复杂度和稳定性

- 18 -

现代信号分析

为代价。

2. 胎儿心电监护

胎儿心电图(Fetal Electrocardiogram，FECG)是研究胎儿心脏电生理活动的一项客观指标，反映了胎儿在孕期中的成长和健康状况。胎儿心电一般是从母亲腹部采集获得，不可避免地会含有许多噪声。首先，加在母体上的工频干扰；其次，母体心电（Maternal Electrocardiogram，MECG），而且与胎儿心电相比，它会强得多；最后，还有夹杂的除心电以外的其他生物电信号。虽然采集获得的信号含有许多种类的噪声，但经前期处理均可滤掉，或者得到很好的抑制，最后剩下的主要是母体心电。因此，我们采用自适应滤波消噪算法把混合信号中的母体心电信号去除。

自适应噪声抵消滤波器的结构框图如图3-2所示：

图4-2 自适应噪声抵消滤波器

上图中抵消器的“信号源”为s+n0，其中s为有用信号，n0为一个与信号s不相关的噪声，抵消器的“噪声源”为噪声n1，n1与信号s不相关，却以某种未知的方式与噪声n0相关，由图可以看出噪声n1经自适应滤波器输出y，再从原始输入s+n0中减去该输出，产生了系统的输出e=s+n0-y。本题中，信号源为母体腹部心电信号，噪声源为母体胸部心电信号MECG，系统输出e为胎儿心电信号FECG。

数据是从MIT-BIH数据库下载，进入网址http://www.physionet.org/cgi-bin/atm/ATM，在左侧Input选项里选择Database为胎儿心电图选项，Record项任选一个，Signals项选择为all，右侧Toolbox选择为Export signals as .mat。如图4.3所示。

图4-3 数据下载

设置完上面步骤后，会找到网页里有一个下载项，ecgca323_edfm.mat，点击下载。这就是对我们有用的数据，MATLAB打开后为一个6*10000的数据，我们只要对第一行和第三行的数据进行处理，分别为母亲胸部和腹部的数据信号。我们分别用LMS算法和RLS算法得到胎儿心电。

胎儿心电监护原理：监测胎儿心电，需要将母亲的心音及背景干扰去除。母亲的心音强度通常是胎儿心音的2倍到10倍，肌肉活动及胎儿运动产生的背景干扰也常常大于胎儿心音的强度。采用自适应对消的方法来监护胎儿心电图，将母亲胸部的信号作为参考输入，主要包括母亲心音和背景噪声，

- 19 -

也有胎儿的心音信号，如图4.4所示（数据为MIT-BIH数据库下载）；从母亲腹部取出的信号作为原始输入，包括母亲心音、胎儿心音和背景干扰，如图4.5所示（数据为下载数据），可以看出，这个系统是一个有信号分量泄露到参考输入端的自适应对消系统。我们希望的结果是对消后母亲的心电信号的强度小于胎儿的心电信号强度，且胎儿的心电信号要比较清楚。 程序如下：

load('ecgca102_edfm.mat'); E=val(1,:)./10000; E1=val(3,:)./10000; thorax=E(1:3000); %胸部信号，参考信号 dn=thorax.'; abdomen=E1(1:3000);% 腹部信号，输入信号 xn=abdomen.'; figure(1); plot((1:3000),dn); xlabel('采样点数');ylabel('幅度'); title('母亲心电图'); figure(2); plot((1:3000),xn); xlabel('采样点数');ylabel('幅度'); title('胎儿心电图'); 母亲胎儿心电图：

图4-4 母亲胸部的心电图，作为期望信号 图4-5 胎儿的心电图，作为输入信号

3.1 用LMS算法得到胎儿心电图

比较图4.4和图4.5可以发现，胎儿的心电信号紧随着母亲的心电信号，这里采样点数选取了3000个，阶数采用20阶（即权系数的个数），步长因子?为改变的参数，分别设置为0.005，0.05，0.5。另外由于设置的阶数为20阶，所以权矢量就有20个，如果查看每一个与迭代次数的关系曲线，就要占很大篇幅，由于权系数（滤波器系数）更新的方程是一样的，所以它们的规律是相似的，这里只显示前三个权系数与迭代次数的关系曲线图。其用LMS算法得到胎儿心电图的MATLAB源代码如下（此处

- 20 -

现代信号分析

以步长为0.005为例）：

load('ecgca102_edfm.mat'); E=val(1,:)./10000; E1=val(3,:)./10000; thorax=E(1:3000); %胸部信号，参考信号 dn=thorax.'; abdomen=E1(1:3000);% 腹部信号，输入信号 xn=abdomen.'; M=20; %阶数 itr=3000;% 迭代次数 en=zeros(itr,1); %误差序列 W=zeros(M,itr);% 权系数 mu=0.5;%步长因子 s=zeros(itr,1);%均方误差 %迭代计算 for k=M:itr x=xn(k:-1:k-M+1); %滤波器M个抽头的输入 y=W(:,k-1).'*x; %滤波器输出 en(k)=dn(k)-y; %第k次迭代的误差 W(:,k)=W(:,k-1)+2*mu*en(k)*x; %权系数更新 s(k)=s(k-1)+abs(en(k)).^2; c(k)=sqrt(s(k)/k);%均方误差 end %求最优时滤波器的输出序列 yn=inf*ones(size(xn)); for k=M:length(xn) x=xn(k:-1:k-M+1); yn(k)=W(:,end).'*x; end figure; subplot(311),plot((1:3000),xn-yn,'r',(1:3000),dn,'b'); legend('胎儿心电图','母亲心电图'); xlabel('采样点数');ylabel('幅度'); subplot(312);plot(c); xlabel('迭代次数');ylabel('均方误差'); subplot(313);plot(W(1,:),'r');hold on; plot(W(2,:));plot(W(3,:),'k'); xlabel('迭代次数');ylabel('滤波器系数'); legend('权系数1','权系数2','权系数3'); 修改步长参数，即mu，分别取0.005，0.05，0.5。运行程序，得到的仿真图形分别如图4-6、图4-7、图4-8所示。

图4-6 步长为0.005时的LMS算法 图4-7 步长为0.05时的LMS算法

- 21 -

图4-8 步长为0.5时的LMS算法

比较上面3幅图，可以很明显的发现，随着步长因子的增加，LMS算法的收敛速度，失调量也随着变化，当步长因子?增大，系统的收敛速度变快，但是相应的失调量也越来越大（如图4.7所示），当步长因子?大到某个数值时，系统不再收敛，也就是系统不再稳定（如图4.8所示）；但是这并不意味着?越小越好，因为步长因子?越小，则系统收敛速度变慢，这在实际应用中也是不容许的。因此参数?的选择应从整个系统要求出发,在满足精度要求的前提下,尽量减少自适应时间。

3.2 用RLS算法得到胎儿心电图

递归最小二乘（RLS）算法，其实现步骤有点类似卡尔曼滤波器，也是通过几个递归公式来循环。采样点数3000，滤波器阶数设置为80，初始化条件设置为：R?1(0)??I,?是小的正数，取0.001。

w(0)?wI?0。选择改变的参数为遗忘因子?，分别取0.97、0.7、0.5，观察滤波效果。同样，由

于滤波器系数过多，这里只观测前三个滤波器系数与迭代次数的曲线。

RLS算法得到胎儿心电图的MATLAB源代码如下（此处以遗忘因子?=0.97为例）： %迭代计算 for k=M:itr x=xn(k:-1:k-M+1);%滤波器M个抽头的输入 K(:,k)=(Rxx*x)./(la+x.'*Rxx*x);%更新滤波器增益 e(k)=dn(k)-W(:,k-1).'*x;%误差序列 W(:,k)=W(:,k-1)+K(:,k)*e(k);%更新权系数 Rxx=(Rxx-K(:,k)*x.'*Rxx)/la;%更新自相关矩阵逆更新 y(k)=W(:,k).'*x;%滤波器输出 e1(k)=dn(k)-y(k); s(k)=s(k-1)+abs(e1(k)).^2; c(k)=sqrt(s(k)/k);%均方误差 end figure; subplot(311), plot((1:3000),xn-y,'r',(1:3000),dn,'b'); legend('胎儿心电图','母亲心电图'); xlabel('采样点数');ylabel('幅度'); subplot(312);plot(c); xlabel('迭代次数');ylabel('均方误差'); subplot(313);plot(W(1,:),'r');hold on; plot(W(2,:));plot(W(3,:),'k'); xlabel('迭代次数');ylabel('滤波器系数'); legend('权系数1','权系数2','权系数3');

- 22 -

现代信号分析

仿真波形分别如图4-9、图4-10、图4-11所示。

图4-9 遗忘因子为0.97时的RLS算法 图4-10 遗忘因子为0.7时的RLS算法

图4-11 遗忘因子为0.5时的RLS算法

比较上面3幅图，根据均方误差与迭代次数的关系曲线图可以发现当?=0.9时，系统在1800点处收敛；而当?=0.5时，系统在500点就收敛了，所以得出结论：遗忘因子?越小，系统的跟踪能力越强，但同时可以发现其对噪声也越敏感；相反，遗忘因子?值越大，系统跟踪能力减弱，但对噪声不敏感，收敛时估计误差也越小。所以，在用普通递推RLS算法时，一定要对?有个准确的取值，才能保证系统性能的最佳状态，所以一般文献要求?的取值范围为0。 .95???0.995

五、维纳滤波器

题目：观测信号为y(k)?s(k)?e(k),其中有信号是s(k)为恒量平稳序列，其统计特征已求得为

?E[s(k)]?0?E[s2(k)]??2s?

噪声e(k)是零均值白噪声，且与有用信号不相关，即

- 23 -

?E[e(k)]?02 ??E[e(i)e(j)]??e?(i?j)

??E[e(k)s(k)]?0求维纳滤波器？（10分）

维纳滤波器实际上就是在最小均方误差条件下探确定滤波器的冲击响应h(n)或系统函数H(z)，也就是求解维纳-霍夫方程的问题。因噪声e(k)是零均值白噪声，且与有用信号不相关,根据已知条件

?E[s(k)]?0?E[s2(k)]??2

s??E[e(k)]?0?2 及 ?E[e(i)e(j)]??e?(i?j)

??E[e(k)s(k)]?0 可求得

rsy(n)?E?s(k)y?(k)?

?E?s(k)s?(k?n)??E?s(k)e?(k?n)??rs(n)??s2ry(n)?E?y(k?n)y?(k)?（s(k)为恒量平稳信号）

?E?s(k?n)?e?(k?n)??E?s(k)?e?(k)??rs(n)?re(n)??s2??e2?(i?j) 则维纳霍夫方程为?Rs?Re?h?rd。

以一阶维纳滤波器为例，则维纳霍夫方程可写为

??s2??e2?s2??h(0)???s2???2? ??222???s??e??h(1)???s???s 求解该方程得

?h(0)?1??h(1)?2?2??2??se 维纳滤波器为：

??s2??2? ??s?H(z)? 均方误差

?s22???2s2e(1?z?1)

- 24 -

现代信号分析

?(k) ??Es(k)?s?2??s2?e2 ?rs(k)??h(k)r(k)?222?s??ek?01?sy

六、FIR维纳滤波器

题目：假定下图所需的信号的

d(n)是一个正弦序列

d(n)=

sin(n?0??),?0在0~?之间任取一值，噪声序列?1?n?和?2?n?都是AR(1)过程，分别由如下

的一阶差分方程产生：

v1(n)??v1(n?1)?e(n)v2(n)??0.6v1(n?1)?e(n)

式中0???1,e(n)是零均值，单位方差的白噪声，与d(n)不相关，取不同的?值，如0.1，0.9，试采,2个不同阶数的FIR维纳滤波器对信号x(n)?d(n)??1(n)进行滤波。并分析其结果。(12分)

维纳滤波器（Wiener filter）是由数学家维纳（Rorbert Wiener）提出的一种以最小平方为最优准则的线性滤波器。其数字系统描述如下图所示，它是要从噪声观测x(n)=d(n)+v(n)中恢复信号d(n)。

d(n)δ(n)W(z)^d(n)+-+e(n)

图6-1 维纳滤波器的一般框图

题目中所示系统的Wiener-Hopf方程可以推到如下：

将v2(n)作为维纳滤波器的输入，以估计噪声v1(n)，则Wiener-Hopf方程为：其中Rv2Rv2??rv1v2，是v2(n)的自相关矩阵，rv1v2是期望信号v1(n)和滤波器输入v2(n)之间的互相关矢量。对互相关量，可推导

- 25 -

rv1v2 =E{ v1(n) v2*(n-k)}= E{ [x(n)-d(n)] v2*(n-k)} = E{ x(n) v2*(n-k)}- E{ d(n) v2*(n-k)}

题目中v2(n)与d(n)不相关，则上式第二项为零，上式简化为

rv1v2 = E{ x(n) v2*(n-k)}

因此，Wiener-Hopf方程变为：

Rv2 ω = rxv2

式中Rv2是v2(n)的自相关矩阵，ω是维纳滤波器系数矢量，rxv2是x(n)和v2(n)的互相关矢量。因此，要求解出维纳滤波器，需要先求出维纳滤波器系数矢量ω。

根据上述分析，在matlab中，首先利用白噪声产生v1(n)和v2（n）序列，之后求解出v2(n)的自相关矩阵和x(n)和v2(n)的互相关矢量，然后求出维纳滤波器系数矢量ω。再利用已求解的维纳滤波器进行题目所示的噪声滤除工作。取ω0=0.05，具体求解过程如下程序所示：

clear; clc; N=512; n=1:1:N; en=randn(1,length(n)); %白噪声 d=sin(0.05*n*pi+pi/3); k=2; %维纳滤波器的阶数 v1=filter([1],[1 -0.1],en); v2=filter([1],[1 0.6],en); x=d+v1; Rv2=xcorr(v2,'unbiased'); %计算序列的自相关函数 Rv1v2=xcorr(v1,v2,'unbiased'); %计算序列的互相关函数 for i=1:k; %构造自相关矩阵和互相关矩阵, Rv21(i,:)=Rv2(N-i+1:N-i+k); Rv1v21(i,:)=Rv1v2(N+i-1); end; h=inv(Rv21)*Rv1v21; v11=filter(h',1,v2); v1_1=filter(h',1,v2); e=x-v11; subplot(4,1,1),plot(n,d,'b'),title('输入信号波形') subplot(4,1,2),plot(n,x,'b'),title('信号加白噪声') subplot(4,1,3),plot(n,v1_1,'r'),title('白噪声估计') subplot(4,1,4),plot(n,e,'g'),title('恢复后信号波形') e1=e-d; a=sum(e1.^2) 当?- 26 -

=0.1，维纳滤波器阶数为2阶时，滤波结果如下图6-2所示。

现代信号分析

图6-2

?=0.1， 2阶维纳滤波器滤波效果

?=0.1，维纳滤波器阶数为5阶时，滤波结果如下图6-3所示。

图6-3?=0.1， 5阶维纳滤波器滤波效果

?=0.9，维纳滤波器阶数为20阶时，滤波结果如下图6-4所示。

图6-4?=0.9， 20阶维纳滤波器滤波效果

?=0.9，维纳滤波器阶数为50阶时，滤波结果如下图6-5所示。

- 27 -

当

当当

图6-5

?=0.9， 50阶维纳滤波器滤波效果

根据上述分析和结果图进行对比可以得出，当?=0.1时，维纳滤波器阶数为5时，维纳滤波器滤波得到的信号幅值与相位与原信号就有较好的一致性。但当当?=0.9时，维纳滤波器阶数为20时，维纳滤波器滤波得到的信号幅值与相位与原信号都没有太好的一致性，只有阶数为50时，滤波得到的信号幅值与相位与原信号才有较好的一致性。可以看出维纳滤波器阶数越高，滤波效果更好一些，在选择50阶滤波时，曲线平滑且接近原波形曲线，但实际情况下，并不是阶数越高越好，当滤波阶数过大时，维纳滤波器在仿真计算过程中会引入舍入误差，当舍入误差逐渐积累，会导致高计算量的高阶情况下均方误差增大，引起滤波不理想。但当超过一定阶数时，阶数的增加对滤波效果的改善已经基本可以忽略，此时单纯靠提高阶数的办法不能改进滤波效果了。说明此时维纳滤波器的效能已经达到最佳，若要再改善效果，必须选用结构更先进的滤波器。

七、AR/MA滤波器

题目：用雷达测量地球和月球之间的距离d，测量过程用下列方程描述

x(n)?d?w(n)

其中w(n)是均值为零，方差为?w的白噪声序列，它表示测量误差。为了提高测量精度，现采用以下两种滤波器分别对进行处理x(n)，试比较其方差?w的大小。 滤波器１： 滤波器２：

式中0?a?1。（10分）

y(n)?ay(n?1)?(1?a)x(n) y(n)?ax(n)?(1?a)x(n?1)

222 由于测量过程x(n)为测量真值d与零均值、方差为?w的白噪声序列，因此，在滤波前，x(n)的均值为d，方差为白噪声序列方差?w，即：

- 28 -

2现代信号分析

E[x(n)]?d VAR[x(n)]??w2本题是对测量过程x(n)经过两种滤波器后信号序列的方差的比较。题中所给出的滤波器1和滤波器2，分别为AR和MA模型，两个滤波器中a值未确定，但给出了相应范围，即0?a?1，为了对这个范围内滤波器1和滤波器2的方差进行比较，在MATLAB程序中，需要设定一个参数，其范围为0至1，采样频率为1000，并将其存入一个数组中，以比较不同a值下的方差大小，即两种滤波器对原信号的滤波效果。

另外为了得出白噪声序列方差对滤波后方差的影响，可以通过输入白噪声不同的输入值获得，本题取方差为0.04观测其效果，Matlab程序见如下：

clear;clc; N=1000; var0=input('输入白噪声方差：'); init=2055615866; randn('seed',init); w=sqrt(var0)*randn(1,N); d=3800;%距离 x=d+w;%输入信号 y1(1)=d; y2(1)=d; for i=1:N a=i/N; for j=2:N y1(j)=a*y1(j-1)+(1-a)*x(j); y2(j)=a*x(j)+(1-a)*x(j-1); end v1(i)=var(x); %输入信号方差 v2(i)=var(y1); %滤波器1后的方差 v3(i)=var(y2); %滤波器2后的方差 end t=(1:N)/N; plot(t,v1,'r',t,v2,'b',t,v3,'g'); xlabel('a值');ylabel('方差'); legend('未滤波','滤波器1','滤波器2'); 运行结果如下图所示： - 29 -

图7-1

?2=0.04时 不同a值的方差比较

如图7-1可知，运用滤波器滤波后的信号方差小于滤波前，即提高测量精度；滤波器1（AR）的信号方差在a值在小于0.2时信号方差和滤波器2（MA）得到的大致相同；但当a值慢慢增大时，大概到点0.53时，滤波器2的信号方差也突然随之增大，而滤波器1（AR）的方差值却始终在减小，因此可知滤波器1的性能要比滤波器2好。因此在选择用滤波器2（MA）时应该选择a值尽量小于其方差开始增大的点，即0.5左右。

下面给出a值为0.3和0.9的情况下，对白噪声序列方差不同时对滤波效果进行了比较，具体程序Matlab程序如下：

clear;clc; N=1000; for i=1:N var0=i/N; init=2055615866; randn('seed',init); w=sqrt(var0)*randn(1,N); d=3800;%距离 x=d+w;%输入信号 y1(1)=d; y2(1)=d; K=900; a=K/N; %确定a值 for j=2:N y1(j)=a*y1(j-1)+(1-a)*x(j); y2(j)=a*x(j)+(1-a)*x(j-1); end v1(i)=var(x); %输入信号方差 v2(i)=var(y1); %滤波器1后的方差 v3(i)=var(y2); %滤波器2后的方差 end t=(1:N)/N; plot(t,v1,'r',t,v2,'b',t,v3,'g'); xlabel('白噪声方差');ylabel('方差'); - 30 -

现代信号分析

legend('未滤波','滤波器1','滤波器2 运行结果见图7-2。

（a）a=0.3 （b）a=0.9

图7-2 不同a值白噪声方差对滤波结果的影响

由图可以看出，两种滤波器下信号方差的变化趋势都为线性增大的状态，当a=0.9时，滤波器1（AR

）性能优于a=0.3， 且明显好于a=0.9时滤波器2（MA），在这种情况下选择滤波器时优先考虑 滤波器1（AR）。

八、卡尔曼滤波器

题目：简述卡尔曼滤波器。并计算下面的例子：（15分）

一连续平稳的随机信号x(t)，自相关函数rx(?)?e??，信号x(t)为加性噪声所干扰，噪声是白噪声，测量值的离散值z(k)为已知，Ts=0.01s,

-3.2,-0.8,-14,-16,-17,-18,-3.3,-2.4,-18,-0.3,-0.4,-0.8,-19,-2.0,-1.2,-11,-14,-0.9,0.8,10,0.2,0.5,-0.5,2.4,-0.5,0.5,-13,0.5,10,-12,0.5,-0.6,-15,-0.7,15,0.5,-0.7,-2.0,-19,-17,-11,-14.

利用卡尔曼滤波，估计信号x(t)的波形。并给出不同的量测噪声对卡尔曼滤波的影响，同时给出迭代次数与均方误差的关系曲线。

1．卡尔曼滤波器简介

卡尔曼滤波是一种高效率的递归滤波器(自回归滤波器), 它能够从一系列的不完全包含噪声的测量中，估计动态系统的状态。

简单来说，卡尔曼滤波器是一个“optimal recursive data processing algorithm（最优化自回归数据处理算法）”。对于解决很大部分的问题，他是最优，效率最高甚至是最有用的。他的广泛应用已经超过30年，包括机器人导航，控制，传感器数据融合甚至在军事方面的雷达系统以及导弹追踪等等。近年来更被应用于计算机图像处理，例如头脸识别，图像分割，图像边缘检测等等。

- 31 -

首先，我们先要引入一个离散控制过程的系统。该系统可用一个线性随机微分方程（Linear Stochastic Difference equation）来描述：

X(k)?A?X(k?1)?B?U(k)?W(k) （8-1） 再加上系统的测量值：

Z(k)?H?X(k)?V(k) （8-2）

U(k)是k时刻对系统的控制量。 上面两式子中，X(k)是k时刻的系统状态，A和B是系统参数，

对于多模型系统，他们为矩阵。Z(k)是k时刻的测量值，H是测量系统的参数，对于多测量系统，H为矩阵。W(k)和V(k)分别表示过程和测量的噪声。他们被假设成高斯白噪声(White Gaussian Noise)，他们的协方差矩阵分别是Q，R（这里我们假设他们不随系统状态变化而变化）。

对于满足上面的条件(线性随机微分系统，过程和测量都是高斯白噪声)，卡尔曼滤波器是最优的信息处理器。下面我们来用他们结合他们的协方差矩阵来估算系统的最优化输出。

首先我们要利用系统的过程模型，来预测下一状态的系统。假设现在的系统状态是k，根据系统的模型，可以基于系统的上一状态而预测出现在状态：

X(kk?1)?A?X(k?1k?1)?B?U(k) (8-3)

U(k) 式(8-3)中，X(kk?1)是利用上一状态预测的结果，X(k?1k?1)是上一状态最优的结果，

为现在状态的控制量，如果没有控制量，它可以为0。

到现在为止，我们的系统结果已经更新了，可是，对应于X(kk?1)的协方差矩阵还没更新。我们用P表示协方差矩阵：

P(kk?1)?A?P(k?1k?1)?A'?Q (8-4)

式(5-4)中，P(kk?1)是X(kk?1)对应的协方差矩阵，P(k?1k?1)是X(k?1k?1)对应的协方差矩阵，A'表示A的转置矩阵，Q是系统过程的协方差矩阵。式子8-1，8-2就是卡尔曼滤波器5个公式当中的前两个，也就是对系统的预测。

现在我们有了现在状态的预测结果，然后我们再收集现在状态的测量值。结合预测值和测量值，我们可以得到现在状态X(k)的最优化估算值 X(kk)：

X(kk)?X(kk?1)?Kg(k)?(Z(k)?H?X(kk?1)) (8-5)

- 32 -

现代信号分析

其中Kg为卡尔曼增益(Kalman Gain)：

Kg(k)?P(kk?1)?H'/(H?P(kk?1)?H'?R) (8-6)

到现在为止，我们已经得到了k状态下最优的估算值X(kk)。但是为了要另卡尔曼滤波器不断的运行下去直到系统过程结束，我们还要更新k状态下X(kk)的协方差矩阵：

P(kk)?(I?Kg(k)?H)?P(kk?1) (8-7)

其中I为1的矩阵，对于单模型单测量，I?1。 2.利用卡尔曼滤波，估计信号的波形 （1）求过程的功率谱密度

将积分按＋?和－?分成两部分进行，信号的自相关函数是rx(?)?e??，则其功率谱为

S0??X(?)??e?e?jd????e??e?j????0d?

e(??j?)?0(??j?)????e?1?j????(1?j?) 0???11?? ?1?j??1?j???21??2

（2） 用复频率s来表示功率谱密度 （???js）

Sx(s)?21??2?11?2??s2s?1?1?s （3）谱分解与Z变换根据冲激响应不变原理，将连续系统变成离散系统，即

Sx(z)?m??1?emTz?m?????e?mTz?me?Tzz1?e?2T?m?01?e?Tz?z?e?T?(1?e?Tz?1)(1?e?Tz) 其中z?ejwT，T?TS?0.01

S2?1(1?e?2T)z?1根据等式zz?1x(z)??wB(z)B(z)?(1?e?Tz?1)(1?e?Tz)得 B(z)?1?e?Tz?1

变换到时域得：x(n)?e?Tx(n?1)?w(n?1),

因此A(k)?e?T?e?0.01?0.99，

- 33 -

所以得出状态方程x(k?1)?0.99x(k)?w(k)，即B=1， 又因为量测方程y(k)?x(k)?v(k)，所以C(k)?1，D(k)=0，

2T=0.01s，Q??w?1?e?2T?0.02。

（4）下面给出的matlab程序，包含自编卡尔曼滤波算法和matlab自带卡尔曼函数。

在用matlab自带卡尔曼函数前，先调用sys=ss（A,B,C,D,Ts）函数，这里采样周期Ts=0.02s，

A,B,C,D均已求出，得到离散卡尔曼状态模型，然后调用函数[kalmf,L,P,H,E] = kalman(sys,Q,R),设

计离散卡尔曼滤波器。H和E表示系统稳态的修正矩阵和均方误差，把它们代入估计值的更新公式（如

下）就可以根据观测信号得到卡尔曼滤波的估计值了。

S(k)?A*S(k?1)?H*[Y(k)?C*A*S(k?1)]

其中，S(k)是卡尔曼滤波估计值，Y(k)是已知的测量值。

clear all; Y = [-3.2,-0.8,-14,-16,-17,-18,-3.3,-2.4,-18,-0.3,-0.4,-0.8,-19, -2.0,-1.2,-11,-14,-0.9,0.8,10,0.2,0.5,-0.5,2.4,-0.5,0.5,-13,0.5, 10,-12,0.5,-0.6,-15,-0.7,15,0.5,-0.7,-2.0,-19,-17,-11,-14]; N= length(Y); % ----卡尔曼函数 ---------------------------------------- A=0.99;B=1;C=1;D=0; T=0.01; sys = ss(A,B,C,D,T); % 求出离散卡尔曼状态模型 v=std(randn(1,N)); R=v^2; Q=0.02; [kalmf,L,P,H,E] = kalman(sys,Q,R); % H：修正矩阵 E:最小均方误差 S=zeros(1,N); for k=2:N S(k)=A*S(k-1)+H*(Y(k)-C*A*S(k-1)); % 卡尔曼滤波后的估计值 end k=1:1:42; plot(k*0.02,S,'b-',k*0.02,Y,'r'); legend('估计值','测量值'); ylabel('信号幅值'); title('卡尔曼函数滤波结果'); grid on; 运行结果如图8-1所示： - 34 -

现代信号分析

图8-1 卡尔曼滤波估计图

上图是直接用卡尔曼滤波函数进行滤波的效果。在Matlab当中，可以根据状态方程和观测方程，调用dlqe函数，直接求解其最佳递推估计方程，可以利用该方程对输入信号进行滤波，比较两种方式的结果，进行结果验证。结果验证的Matlab程序如下：

% 采用dlqe函数求解kalman滤波器的稳态参数 a=0.99;G=1;Qk=0.02;Rk=1;c=1; [L,P1,P]=dlqe(a,G,c,Qk,Rk); %%%%%%%%% 根据kalman滤波器的稳态参数设计的滤波器 for i=1:size_x(1) xw(i)=X(i); end for i=2:size_x(1) xw(i)=0.99*xw(i-1)+L*(X(i)-0.99*xw(i-1)); end % 结果显示 subplot(312);plot(k*0.02,xk_s(k),'b-',k*0.02,xw(k),'r'); legend('Kalman','Wiener'); xlabel('时间Ts');ylabel('信号幅值');title('Kalman和稳态滤波比较');grid on; 本题中，信号x(t)为加性噪声所干扰，噪声是白噪声，但是白噪声的方差没有给出， 以??0.05为例，Matlab程序得到的结果如图8-2所示。

2

图8-2 卡尔曼滤波结果和稳态滤波结果比较

- 35 -

稳态滤波结果是采用dlqe函数编写的滤波函数和kalman，从图可以看出，卡尔曼滤波在0.3s以前和稳态滤波有一定的差异，但是在0.3s之后，完成学习的过程，基本和采用dlqe函数编写的滤波函数重合，本题编写的kalman滤波函数可以完成滤波的工作，结果正确。由于白噪声的方差没有给出，下面比较不同白噪声方差（?=0.05、0.5、1）下的结果，结果图如图8-3所示。

2

图8-3 不同?下测量信号、kalman滤波信号 稳态滤波信号比较以及测量信号和kalman滤波信号的偏差

2

从图8-3中，可以看到，随着观测方程的方差?的变大，对于原始波形尖峰的消减越大，得到的曲线也越光滑，kalman滤波器的学习时间也越长。

2

图8-4 不同方差下迭代次数与均方差曲线

- 36 -

现代信号分析

图8-4给出了方差分别为0.05、0.5、1时的迭代次数与均方差曲线图。由图可知当量测噪声方差越小时，迭代越快，所得的均方差也越小。结合图8-3可知量测噪声方差越小的话，估计出来的信号量也越接近实际值。

九、小波分析方法及应用

题目：简要介绍小波分析方法及消噪中的应用，并举一例加以说明（建议采用实验数

据，对阈值的选择要给出依据）。（10分）

1. 小波分析原理

小波分析方法是一种窗口大小（即窗口面积）固定但其形状可改变，时间窗和频率窗都可改变的时频局部化分析方法。即在低频部分具有较高的频率分辨率和较低的时间分辨率，在高频部分具有较高的时间分辨率和较低的频率分辨率，所以被誉为数学显微镜。正是这种特性，使小波变换具有对信号的自适应性。原则上讲，传统上使用傅里叶分析的地方，都可以用小波分析取代。小波分析优于傅里叶变换的地方是，它在时域和频域同时具有良好的局部化性质。

小波函数的定义：设??t?为平方可积函数，即??t??L(R)，若其傅里叶变换????（????是

2????t?的傅里叶变换）满足

?(?)C???R?2?d???

称??t?为一个基本小波或母小波（Mother Wavelet），并称上式为小波函数的允许条件。与标准的傅立叶变换相比，小波分析中用到的小波函数不具有唯一性，对于一个时频分析问题，如何选者最佳的小波基函数是一个重要的问题。常用的小波函数有Haar小波、dbN小波、Morl小波、Mexh小波、Meyer小波等，不同的小波函数对应不同的尺度函数和性能。从下图中可以看出小波变换与傅立叶变换在时频窗口特性上有很大的不同，更显示了上述小波变换的特点。

?t 图9-1 小波变换的时频分析窗

小波变换的多分辨率分析实际上就是对一个频带信号进行低频分解，对每一步分解出来的低频部分在分解，使频率分辨率越来越高，其目的是构造一个理想的正交小波基。小波包分析实际上就是对

- 37 -

与多分辨率分析没有分解的高频信号也进行逐层分解，进一步提高时频分辨率。小波分析地这些原理与特点与测控领域中的滤波原理非常相似，常常被用于信号噪声的消除。

运用小波分析进行一维信号消噪和识别信号中的发展趋势是小波分析的一个重要应用之一。在实际的工程应用中，所分析的信号可能包含噪声和一些不应有的趋势项，对这种信号，首先需要作信号的预处理，将信号的噪声和趋势项去处，提取有用信号。

2. 小波分析的应用

在实际工程应用中，所分析的信号可能包含许多尖峰或突变部分，并且噪声也不是平稳的白噪声，对这种信号进行分析，进行去噪处理，传统的傅立叶分析显得无能为力，因为它不能给出信号在某个时间点上的信号变换情况，使得信号在时间轴上的任何一个突变都会影响信号的整个谱图。而小波分析由于能同时在时域，频域中对信号进行多分辨分析，所以能有效地区分信号中的突变部分和信号噪声，从而实现信号的去噪。下面以小波分析在信号去噪方面的应用为例，分析不同小波函数、分解层数、阈值三个方面分析滤波的效果，介绍小波分析在信号处理方面的的应用。处理的数据为热电偶的输出。

对于小波消噪的算法主要分为三步：

(1)用小波变换或小波包变换将含噪数据变换到小波域； (2)确定一个阈值标准并用作用到小波域； (3)重构数据，得到消噪后的信号；

（1）不同小波函数

小波函数种类繁多，不同的小波函数适用不同的情况，下面以常用的五种小波函数Daubechies小

波、Coiflets小波、Symlets小波、正交小波、双正交小波对本实验数据进行处理，选择最合适的小波函数。

clear all; fid = fopen('热电偶数据.txt'); A = fscanf(fid,'%f',[1 inf]); meg_sensor=A; meg_sensor=meg_sensor(1:5000);lengtha=size(meg_sensor); %%%%%%%%%%%%%%%%%%%%%%%% 小波分解 [c_db,l_db]=wavedec(meg_sensor,6,'db5'); [c_coif,l_coif]=wavedec(meg_sensor,6,'coif5'); [c_sym,l_sym]=wavedec(meg_sensor,6,'sym5'); [c_bior,l_bior]=wavedec(meg_sensor,6,'bior5.5'); [c_rbio,l_rbio]=wavedec(meg_sensor,6,'rbio5.5'); %%%%%%%%%%%%%%%%%%%%%%% 强制消噪 meg_sensor_db=wrcoef('a',c_db,l_db,'db5',6); - 38 -

现代信号分析

meg_sensor_coif=wrcoef('a',c_coif,l_coif,'coif5',6); meg_sensor_sym=wrcoef('a',c_sym,l_sym,'sym5',6); meg_sensor_bior=wrcoef('a',c_bior,l_bior,'bior5.5',6); meg_sensor_rbio=wrcoef('a',c_bior,l_bior,'rbio5.5',6); %%%%%%%%%%%%%%%%%%%%%% 结果显示 %%%%%%%%%%%%%%%%%%%%%%% 原始信号和强制滤波之后的信号比较 figure(1);k=1:1:lengtha(2); subplot(321); plot(k*0.00025,meg_sensor(k)); xlabel('时间T/s');ylabel('信号幅值');title('原始信号');grid on; subplot(322); plot(k*0.00025,meg_sensor_coif(k)); xlabel('时间T/s');ylabel('信号幅值');title('Coiflets小波滤波之后的波形');grid on; subplot(323); plot(k*0.00025,meg_sensor_sym(k)); xlabel('时间T/s');ylabel('信号幅值');title('Symlets小波滤波之后的波形');grid on; subplot(324); plot(k*0.00025,meg_sensor_bior(k)); xlabel('时间T/s');ylabel('信号幅值');title('Biorthogonal小波滤波之后的波形');grid on; subplot(325); plot(k*0.00025,meg_sensor_rbio(k)); xlabel('时间T/s');ylabel('信号幅值');title('Reverse Biorthogonal滤波之后的波形');grid on; subplot(326); plot(k*0.00025,meg_sensor_db(k)); xlabel('时间T/s');ylabel('信号幅值');title('Daubechies小波滤波之后的波形');grid on;

上述程序得到的结果和原始结果如图6.2所示。

图9-2原始信号以及不同小波滤波的结果曲线

从图中可以看出，五种小波函数的滤波效果都比较好，都可以很好的得到平滑的曲线，综合比较

- 39 -

滤波效果和算法，我们选择Daubechies小波处理的加速度传感器输出信号。

（2）不同分解层数

小波分析当中，对信号进行分解，第一层分解将原始信号分解为高频、低频两部分，第二层分解

将第一层分解的低频信号再分解为高频、低频两部分，依次分解下去。采用Daubechies小波进行不同层的分解，采用强制滤波的方法进行信号重构，所编写的Matlab程序如下所示。

clear all; fid = fopen('热电偶数据.txt'); A = fscanf(fid,'%f',[1 inf]);meg_sensor=A;lengtha=size(A); %%%%%%%%%%%%%%%%%%%%%%%% 小波分解 [c3,l3]=wavedec(meg_sensor,3,'db6'); [c4,l4]=wavedec(meg_sensor,4,'db6'); [c5,l5]=wavedec(meg_sensor,5,'db6'); [c6,l6]=wavedec(meg_sensor,6,'db6'); [c7,l7]=wavedec(meg_sensor,7,'db6'); %%%%%%%%%%%%%%%%%%%%%%%% 强制消噪（仅使用cA5重构信号） meg_sensor_filter3=wrcoef('a',c3,l3,'db6',3); meg_sensor_filter4=wrcoef('a',c4,l4,'db6',4); meg_sensor_filter5=wrcoef('a',c5,l5,'db6',5); meg_sensor_filter6=wrcoef('a',c6,l6,'db6',6); meg_sensor_filter7=wrcoef('a',c7,l7,'db6',7); %%%%%%%%%%%%%%%%%%%%%%% 原始信号和强制滤波之后的信号比较 figure(2);k=1:1:lengtha(2); subplot(321); plot(k*0.00025,meg_sensor(k)); xlabel('时间Ts');ylabel('信号幅值');title('原始信号');grid on; subplot(322); plot(k*0.00025,meg_sensor_filter3(k)); xlabel('时间Ts');ylabel('信号幅值');title('分解三层滤波之后的波形');grid on; subplot(323); plot(k*0.00025,meg_sensor_filter4(k)); xlabel('时间Ts');ylabel('信号幅值');title('分解四层滤波之后的波形');grid on; subplot(324); plot(k*0.00025,meg_sensor_filter5(k)); xlabel('时间Ts');ylabel('信号幅值');title('分解五层滤波之后的波形');grid on; subplot(325); plot(k*0.00025,meg_sensor_filter6(k)); xlabel('时间Ts');ylabel('信号幅值');title('分解六层滤波之后的波形');grid on; subplot(326); plot(k*0.00025,meg_sensor_filter7(k)); xlabel('时间Ts');ylabel('信号幅值');title('分解七层滤波之后的波形');grid on; - 40 -

现代信号分析

图 9-3 采用db6小波进行不同层分解得到的滤波之后的信号

从图中可以看出，原始信号中含有较大的噪声，而经过小波滤波处理之后，随着分解层数的增加，滤波效果越来越好。不过，在滤波层数在三层以上时，滤波效果已经没有明显的变化，因此，三层分解最适合于本数据。 （3）阈值对滤波的影响

采用Daubechies小波中的db6，将信号作3层分解，根据分解之后的信号的波形判断每一层的阈值，最后根据阈值重构信号，得到滤波之后的信号，所编写的matlab程序如下。

clear all; %A=importdata('热电偶数据',' '); fid = fopen('偏30仰30.txt');A = fscanf(fid,'%f',[6 inf]);meg_sensor=A(5,:)'; meg_sensor=meg_sensor(1:5000); %%%取1~5000处的数值进行处理，采样频率4K lengtha=size(meg_sensor); %%%%%%%%%%%%%%%%%%%%%%%% 小波分解 [c,l]=wavedec(meg_sensor,5,'db6'); cA3=appcoef(c,l,'db6',3); %%提取小波分解的低频系数 cD3=detcoef(c,l,3); %%提取第三层的高频分量 cD2=detcoef(c,l,2); %%提取第二层的高频分量 cD1=detcoef(c,l,1); %%提取第一层的高频分量 %%%%%%%%%%%%%%%%%%%%%%%% 画出各个分量的系数图 figure(1); subplot(321);plot(cD1);title('cD1');grid on; subplot(322);plot(cD2);title('cD2');grid on; subplot(323);plot(cD3);title('cD3');grid on; subplot(324);plot(cA3);title('cA3');grid on; %%%%%%%%%%%%%%%%%%%%%%%% 强制消噪（仅使用cA3重构信号） meg_sensor_filter=wrcoef('a',c,l,'db6',3); - 41 -

%%%%%%%%%%%%%%%%%%%%%%% 原始信号和强制滤波之后的信号比较 figure(2);k=1:1:lengtha(1); subplot(221)plot(k*0.00025,meg_sensor(k)); xlabel('时间Ts');ylabel('信号幅值'); title('原始信号');grid on; subplot(222);plot(k*0.00025,meg_sensor_filter(k));xlabel('时间Ts');ylabel('信号幅值'); title('强制滤波之后的信号');grid on; 分解之后各个层的系数图如图9-4所示。

图9-4各个层之间的系数结果

直接采用强制消噪得到的结果如左下图所示。从图中可以看出，滤波得到的信号非常平滑。但是如果原始信号如右图所示的振动信号，在1.6s左右的高频振动信号，也完全被滤掉了，使得滤波之后的信号丢失了部分重要的信号。

图9-5原始信号和强制消噪之后的信号比较

为了得到很好的滤波效果并且不丢失所需要的高频分量，下面分别采用默认阈值和自定义软阈值的方法进行消噪。

- 42 -

现代信号分析

图9.6 原始信号和强制消噪、默认消噪、软阈值消噪信号的比较

从图中可以看出强制消噪、默认消噪、软阈值消噪之后的信号，可以得到，强制消噪丢失了部分信息，直接采用默认阈值进行滤波，可以很好的保留所需的高频信号，但是，在低频部分滤波效果不好，采用自定义软阈值的方式，可以兼顾滤波效果和高频分量两部分，滤波效果较好，但需要对信号不同层的分量进行分析，确定有用部分。

十、其他现代信号分析方法

题目：阅读上课内容以外的其他现代信号方法相关文章，请对其内容进行总结，从技术思路、主要优缺点和您的个人评价三个角度加以阐述。（5分）

结合所阅读的文献，在此将介绍独立分量分析在盲信号处理中的应用。首先简要的介绍下盲信号处理的一些基本知识，并进一步分析独立分量分析在盲信号处理中的应用。

1.盲信号处理

盲源分离是指在信号的理论模型和源信号无法精确获知的情况下,如何从混迭信号(观测信号)中分离出各源信号的过程。盲源分离和盲辨识是盲信号处理的两大类型。盲源分离的目的是求得源信号的最佳估计，盲辨识的目的是求得传输通道混合矩阵。其原理框图如图10-1所示。

- 43 -

10-1 盲信号处理原理图

盲信号处理的主要研究方向有：

（1）盲信源分离（BSS）。源信号、传输通道特性未知，由观测信号和源信号先验知识估计出源信号各个分量过程；

（2）独立分量分析（ICA)。针对独立源信号混合的各分量问题提出独立分量分析。在源信号相互独立的假设下，ICA与BSS具有相同模型；

（3）多通道盲解卷积（BD）/盲均衡（BE）。针对源信号卷积混合的分离问题。混合系统特性参数未知，由观察的混合信号重构源信号，可应用到时变系统及非最小相位系统。 盲信号处理方法和主要思路

（1）HOS：高阶统计量衡量信号的独立性和高斯性，或稀疏性（ICA）；

（2）SOS：有时序结构用二阶统计量（SOS）即可，不能分离具有相同功率谱形状或独立同分布信号；

（3）NS+SOS：利用非平稳信息和SOS结合，能够分开功率谱形状相同的源信号。但若非平稳性也相同就不可以分离；

（4）STF多样：运用信号不同多样性：时域多样性，频域多样性，空域多样性。 TDMA， FDMA， SDMA

2.独立分量分析（ICA）的基本原理

1994年, Comon系统地分析了瞬时盲信号分离问题，同时明确地引入独立分量分析(ICA) 的概念, 证明了只要恢复出混合信号中各个信号之间的相互独立性, 就可以完成对源信号的分离。可以说, Comon的工作实际上使得对盲信号分离算法的研究变成了对独立分量分析的代价函数以及其优化算法的研究, 使得以后的算法设计开始有了明确的理论依据。

独立分量分析的含义是把信号分解成若干个互相独立的成分，它是为了解决盲信号分离的问题而发展起来的。如果信号本来就是由若干独立信源混合而成的，ICA就能恰好把这些信源分解开来。故在一般的文献中通常把ICA等同于BSS，ICA不同于主分量分析把目光投注于信号的二阶统计量，研究信号间的相关关系，而是基于信号的高阶统计量，研究信号间的独立关系。

因此ICA的前提是：假设S中各分量相互独立；零均值，且方差为1。以多导信号处理为基础，必须借助于一组把信源按不同比例组合起来的多通道信号同步观察。多导信号包括：主分量分析（PCA）和奇异值分解（SVD）。

ICA的简单思路ICA的任务明确为：在S，A均未知的情况下，求B，使Y=BX是S的最优逼迫。

图10-2 ICA分离

服从的基本原则: (1)非线性去相关。求B，使任意两输出yi, yj(i≠j)不相关; 且经非线性变换g(yi), h(yi)也不相关（高阶统计量）。 (2) 使输出尽可能非高斯化。Y的非高斯性的每个局部极大值都给了

- 44 -

现代信号分析

一个独立分量。

ICA模型模型在分为两种，无噪声模型和有噪声模型。 （1）无噪声的ICA模型

ICA作为生成模型的估计给定随即变量的一组观测本标号，假设它们有独立成分线性混合而产生：

?x1(t)??s1(t)??????x2(t)??s2(t)??x3(t)??s3(t)?????x(t)s(t)??4?4? (10-1) ?=A?x1(t),x2(t),x3(t)…xn(t)，其中t是时间或者样

式中，A是某个未知矩阵。

用向量-矩阵符号方式表示通常比上面的求和表达式更为方便。用随机向量x来表示混合向量，其元素分别为x1,...,xn，同样地，用s来表示元素s1,...sn，用矩阵A表示那些混合系数aij。所有的向量都理解为列向量；这样x或者称x的转置就是一个行向量。利用向量和矩阵符号表示，混合模型可以写为：

x?As (10-2)

有时我们需要使用矩阵A中的列向量，如果将其表示为,则模型也可以写为：

Tx??aisii?1n (10-3)

（2）有噪声的ICA模型

将基本的ICA模型扩展到有噪声的情形，并且假设噪声是以加性噪声形式存在的。这是一个相当现实的假设，因为加性噪声是因子分析和信号处理中通常研究的标准形式，具有简单的噪声模型表达方式。因此，噪声ICA模型可表示为：

x?As?n (10-4)

式中，

n??n1,...nn?T是噪声向量。

信号源噪声，即直接添加到独立成分（即信号源）上的噪声。信号源噪声可用与式（10-1）稍有差别的下式来表示：

x?A(s?n) (10-5)

实际上，如果可以直接考虑带噪声的独立成分，那么可将此模型写为：

x?As (10-6)

~可以看出，这就是基本的ICA模型，只是独立成分本身变了。

3.独立分量分析的算法实现

独立分量分析的算法在近几年的发展中出现了很多，在此将简单介绍几种

（1）负熵的FastICA算法(FastICA using negentropy algorithm)

负熵源于信息论中熵的概念，定义如下：

- 45 -

本文来源：https://www.bwwdw.com/article/kwfa.html