现代信号处理
现代信号处理
现代信号处理
现代信号处理是对信号进行数字化处理的一种技术，它使用数字信
号处理算法来分析、修复、增强或压缩信号。
现代信号处理技术广
泛应用于通信、音频处理、图像处理、生物医学工程、雷达和声纳
等领域。
现代信号处理的基本步骤包括信号采集（模拟信号转换为数字信号）、滤波、采样、量化和编码。
滤波可以用于去除信号中的噪声
或不需要的成分，采样和量化将连续的信号转换为离散的数据点，
编码则将离散的数据点转换为数字形式，方便存储和传输。
现代信号处理算法包括傅里叶变换、小波变换、自适应滤波、功率
谱估计以及各种滤波器设计方法等。
傅里叶变换可以将信号从时域
转换为频域，从而可以分析信号的频谱特性；小波变换可以将信号
分解成不同的频率分量，实现信号的多分辨率分析；自适应滤波可
以根据信号的特性自动调整滤波器的参数，以适应不同的环境条件。
1
现代信号处理技术在通信领域广泛应用，例如调制解调、信道编码、多址接入等；在音频处理中，可以实现音频降噪、语音识别和语音
合成；在图像处理中，可以实现图像去噪、边缘检测和数字图像压缩；在生物医学工程中，可以实现生物信号的特征提取、滤波和分析；在雷达和声纳中，可以实现目标检测、目标跟踪和图像重建。
总之，现代信号处理技术为信号分析和处理提供了一种高效、准确
和灵活的方法，为我们获取有用的信息、改善信号质量和实现更复
杂的信号处理任务提供了重要的工具。
2。
现代信号处理盲
SCA利用信号的稀疏性进行盲信号处理，通过寻找观测信号中的稀疏 成分来恢复源信号。
非负矩阵分解（NMF）
NMF是一种基于非负性约束的矩阵分解方法，可用于盲信号处理和特 征提取。
深度学习
近年来，深度学习在盲信号处理领域取得了显著进展，通过训练深度 神经网络模型来实现盲信号处理和源信号分离。
01
信号处理基础
信号定义与分类
信号定义
信号是传递信息的物理量，可以 是电信号、光信号、声信号等。 在信号处理中，主要研究电信号 的处理。
信号分类
根据信号的性质和特征，信号可 分为模拟信号和数字信号、连续 时间信号和离散时间信号、确定 性信号和随机信号等。
线性时不变系统
线性系统
线性时不变系统的性质
线性系统是指系统的输出与输入之间满 足线性叠加原理，即输出的总响应等于 各输入单独作用时产生的响应之和。
线性时不变系统具有稳定性、因果性、 可逆性、可交换性等性质，这些性质 在信号处理中具有重要意义。
时不变系统
时不变系统是指系统特性不随时间变 化，即输入信号的时移不会导致输出 信号的时移。
频域分析与变换
REPORT
THANKS
感谢观看
CATALOG
DATE
ANALYSIS
SUMMAR Y
信号失真比（SDR） 反映输出信号相对于原始信号的失真程度，值越 高表示分离效果越佳。
3
源信号与估计信号的相关系数
通过计算源信号与估计信号之间的相关系数，评 估分离算法对源信号的恢复程度。
计算复杂度评估
算法运算量
统计算法在执行过程中所需的乘法、加法等基本运算次数，以评 估其计算复杂度。
算法执行时间
现代信号处理_完美版PPT
•
测量信号v(n)是均值为零，方差为
2 v
的高斯白噪声；
且v(n)与信号x(n)统计无关，即v(n)不影响信号的谱形状
故有
S y ( y ) S x (x ) v 2 u 2 H () 2 v 2 R u ( m y ) E [ u ( n ) y ( n m ) ] u 2 h ( m )
2
高阶谱估计
➢ 研究的必要性 ➢ 高阶统计量 ➢ 高阶谱 ➢ 高阶累积量和多谱的性质 ➢ 三阶相关和双谱及其性质 ➢ 基于高阶谱的相位谱估计 ➢ 基于高阶谱的模型参数估计 ➢ 多谱的应用
参考:《现代数字信号处理》(184-199；204-205)
3
研究高阶谱的必要性
❖ 关于模型参数估计问题
• 所谓模型参数估计，就是根据有限长的数据序列(如模 型输出端所观测到的信号y(n)来估计图中随机信号模型 的参数，)
i1
i1
即不同ARMA过程具有相同形状的功率谱。这一特性 称为相关函数的多重性或模型的多重性。
9
随机信号的高阶特征（续）
两个具有零均值和相同方差的高斯白色噪声和 指数分布白色噪声显然是不同的随机过程，但它 们的功率谱相同。
用这样两个白色噪声激励同一个ARMA模型，产生的 两个ARMA过程显然是不同的随机过程，但它们的
• 与前面所述不同之处在于：这里考虑了观测过程所引 入的噪声v(n).
v(n)
u(n)
H(z)
x(n) ∑
y(n)
(h(n))
4
研究高阶谱的必要性
❖ 基于二阶统计量的模型参数估计方法的缺陷
• 前述模型参数估计方法中，估计得到的模型参数仅与 信号的自相关函数或功率谱包络相匹配；其功率谱不 含信号的相位特性，亦称盲相。即
现代信号处理的几个边沿问题
3. 信号分析方法只限于二阶矩特性和傅氏频谱。
4. 傅里叶变换的困境
○ 在信号分析和故障诊断技术等领域中，以前最为普遍
○ 是利用快速傅里叶变换 (FFT) 的频域分析法，这种方法
MATLAB 仿真见图1 。
图1 正弦波与回 声信号叠加的波 形和时谱形状
衬底1
Signal in time domain 1
0.5
0
-0.5
-1
0
0.5
1
1.5
Time/s
Cepstrum of signal 1
0.5
0
-0.5
-1
0
0.5
1
1.5
Time/s
(2) 功率频谱(不是功率时谱)
短时: 小时间 区间。
衬底1
应用举例: 开关电源 传导干扰信号的短时 分形维数模糊控制滤 波
基于短时分形维数的模糊控制滤波方法， 对开关电源传导干扰信号中的噪声进行滤 波。该方法提出了网络分形维数和短时分 形维数的新算法，并讨论了模糊控制滤波 方法中的模糊控制参数的选取算法。基于 虚拟仪器(VI) LabVIEW 6.i平台上对开关 电源传导干扰信号进行实时检测。经过信 号处理，该系统还具有信噪分离、测量传 导干扰功率谱等功能。结果表明，该方法 滤波效果良好。
Tga,t0a 1 f(t)g t at0 dt
1 g t t0 a a
其中小波 是将具有局部特性的小 波函数g(t)通过平移和尺度变换（放大倍数为1/a）而构成的。参
数a具有时间的量纲，也称 为小波尺度;f(t)为被处理的信号。 小波函数g(t)称为小波母函数,有多种,以便 适应各种非平稳信号的检测。当对信号进行小波 变换时，其局部化特性与所选取小波函数有关， 因此，要根据信号的特征选择适当的小波母函数 才能获得满意的检测效果。
现代信号处理
现代信号处理现代信号处理⼀信号分析基础傅⾥叶变换的不⾜：()()1()()2j t j tX j x t e dtx t X j e d π∞-Ω-∞∞Ω-∞Ω==ΩΩ??1.不具有时间和频率的“定位”功能；2.傅⾥叶变换对于⾮平稳信号的局限性；3.傅⾥叶变换在分辨率上的局限性。
频率不随时间变化的信号，称为时不变信号（⼜称为平稳信号），频率随时间变化的信号称为时变信号（⼜称为⾮平稳信号），傅⾥叶变换反映不出信号频率随时间变化的⾏为，只适合于分析平稳信号。
⽽我们希望知道在哪⼀时刻或哪⼀段时间产⽣了我们所要考虑的频率，现代信号处理主要克服傅⾥叶变换的不⾜，这些⽅法构成了现代信号处理。
分辨率包括频率分辨率和时间分辨率，含义是指对信号能作出辨别的时域或频域的最⼩间隔。
分辨率的好坏⼀是取决于信号的特点，⼆是取决于信号的长度，三是取决于所⽤的算法。
克服傅⾥叶变换不⾜的主要⽅法有：⽅法⼀：STFT （Short Time Fourier Transform ）⽅法⼆：联合时频分析Cohen 分布，联合时频分析Wigner 分布⽅法三：⼩波变换⽅法四：信号的⼦带分解，将信号的频谱均匀或⾮均匀地分解成若⼲部分，每⼀个部分都对应⼀个时间信号。
⽅法五：信号的多分辨率分析，与⽅法四类似，为了适应在不同频段对时域和频域分辨率的不同要求，可以将信号的频谱做⾮均匀分解。
明确概念：时间中⼼、时间宽度、频率中⼼和频带宽度信号能量：2221()()()2E x t x t dt X j d π===ΩΩ<∞??时间中⼼：21()()t t x t dt Eµ=频率中⼼：21()()2x d EµπΩ=ΩΩΩ? 时间宽度：22201()()t t t x t dt E ∞-∞=-频率宽度：22221=()2X d Eπ∞Ω-∞ΩΩΩ-Ω? 时宽和带宽：2,2t T B Ω=?=?品质因数=信号的带宽/信号的频率中⼼。
现代信号处理第八章基于EMD的时频分析方法及其应用
目前EMD方法主要应用于一元信号处理领域，未来研究将拓展其在多元信号处理中的应用，如多 通道信号分析、多维数据融合等。
EMD在复杂系统故障诊断中的应用
复杂系统的故障诊断是信号处理领域的重要研究方向之一，未来研究将探索将EMD方法应用于复 杂系统的故障诊断中，以提高诊断的准确性和可靠性。
01 基于EMD的时频分析方 法概述
EMD方法简介
EMD（Empirical Mode Decomposition）即经验模态分解，是 一种自适应的信号处理方法。
EMD方法能够将复杂信号分解为一系列固有 模态函数（Intrinsic Mode Functions， IMFs），这些IMFs表征了信号在不同时间 尺度上的局部特征。
THANKS FOR WATCHING
感谢您的观看
图像去噪与增强技术
EMD去噪原理
基于经验模态分解（EMD） 的去噪方法通过分解图像信号 为多个固有模态函数（IMF），
有效去除噪声成分。
自适应阈值处理
结合EMD与自适应阈值技术， 实现图像噪声的智能抑制，提
高图像质量。
对比度增强
利用EMD方法对图像进行分 层处理，调整各层对比度，实
现图像整体对比度的增强。
边界效应问题
EMD方法在分解过程中，对信号两端的数据处理存在不确 定性，容易产生边界效应，影响分解结果的精度和可靠性。
发展趋势预测
自适应噪声抑制技术
针对噪声干扰问题，未来研究将更加注重自适应噪声抑制 技术的发展，以提高EMD方法在噪声环境下的性能。
改进EMD算法
为解决模态混叠问题，研究者将致力于改进EMD算法，如引入 掩膜信号、优化筛选过程等，以提高分解的准确性和稳定性。
现代信号处理-胡广书-清华
X ( jΩ)
=
1 2π
<
x(t), e jΩt
>
式中 < x, y > 表示信号 x 和 y 的内积。若 x ， y 都是连续的，则
（1.1.5）
< x, y >= ∫ x(t) y*(t)dt
若 x ， y 均是离散的，则
< x, y >= ∑ x(n) y*(n)
从时域波形还是从频域波形，我们都很难看出该信号的调制类型及其他特点。和图 1.1.1(c)
一样，图 1.1.2(c)也是 x(n) 的时－频分布表示，由该图可明显看出，该信号的频率与时间成
Line ar sca le
Real part
S ignal in time 1
0
-1 |S TF T|2, Lh=48 , Nf=1 92, lin. scale, co ntour, Thld =5%
gt,Ω (τ ) = g(t − τ )e jΩτ
（1.1.8）
来代替傅立叶变换中的基函数 e jΩt ，则
< x(τ ), gt,Ω (τ ) >=< x(τ ), g(t −τ )e jΩτ >
∫= x(τ )g*(t − τ )e− jΩτ dτ = STFTx (t, Ω)
(1.1.9)
该式称为 x(t) 的短时傅立叶变换（Short Time Fourier Transform, STFT）。式中 g(τ ) 是一窗函
愈多。但由傅立叶变换 X ( jΩ) 看不出在什么时刻发生了此种类型的突变。现举两个例子说
明这一概念。 例 1．1．1 设信号 x(n)由三个不同频率的正弦所组成，即
现代信号处理
R x(y)E {x(t)y*(t)}
互协方差函数
C x(y ) E {x ( [ t)x ]y ( [ t )y ] * } Rxy()x*y
互相关系数
xy()
Cxy()
Cxx(0)Cyy(0)
主要性质
1.对零均值随机信号，相关函数与协方差函数
非平稳即不具有广义平稳。 例1.1.1
随机信号的遍历性
均方遍历：一个平稳信号，其n阶矩及较
低阶的所有矩都与时间无关，对所k 有1, ,n
和所有整数 t1,,tk ，恒有
N l i E m 2 N 1 1t N N x (t t1) x (t tk)(t1, ,tk)2 0
及 ，其k阶矩有界，并满足
( t 1 , ,t k ) ( t 1 , ,t k )
广义平稳（协方差平稳、弱平稳）：均值为常 数，二阶矩有界，协方差函数与时间无关。
严格平稳：概率密度函数与时间无关。
3者关系 广义平稳是n=2的n阶平稳； 严格平稳一定广义平稳，反之则不一定；
等价
2. 0 时，自相关函数退化为二阶矩
Rxx(0)E{x(t)2}
3. 0时，协方差函数退化为方差 Cx(x0)Rx(x0)x2
4. R* xx()Rxx() 5. C* xx()Cxx() 6. C x(x)C x(x 0),
R* xy()Ryx()
白噪声
互功率谱密度
定义
P x(yf) Cx(y )ej2fd
互功率谱的实部称为同相谱，虚部称为正交谱。
相干函数
定义 C(f) Pxy(f)
特点
现代信号处理
4.信号的函数表达式为：()()()()sin(2100) 1.5sin(2300)sin(2200)x t t t A t t dn t n t πππ=++++，其中，()A t 为一随时间变化的随机过程，()dn t 为经过390—410Hz 带通滤波器后的高斯白噪声，()n t 为高斯白噪声，采样频率为1kHz ，采样时间为2.048s 。
（1）利用现代信号处理的知识进行信号谱估计；（2）利用现代信号处理知识进行信号的频率提取； （3）分别利用Winner 滤波和Kalman 滤波进行去噪； （4）利用Wigner-Ville 分布分析信号的时频特性。
（1）：利用现代信号处理的知识进行信号谱估计：经典谱估计中两种主要的方法为直接法和间接法，其中间接法则先根据N 个样本数据()x n 的样本自相关函数()()()1*01,01N x n R k x n k x n k M N-==+=⋅⋅⋅∑，，，（4.1）其中1M N ≤<，且()()*x x R k R k -=。
计算样本自相关函数的Fourier 变换，得到功率谱()()Mjk x x k MP R k e ωω-=-=∑（4.2）周期图方法估计的功率谱为有偏估计，可通过加窗来减少其偏差。
定义为 ()()()2101N jn x n P x n c n e NWωω--==∑ （4.3）式中()()122112N n W c n C d NNππωωπ--===∑⎰（4.4）式中，()C ω是窗函数()c n 的Fourier 变换。
功率谱估计程序为： clear clcclose all hidden sf=1000;nfft=2048; t=0:1/1000:2.047;A=normrnd(0,1,1,2048); N=wgn(1,2048,1); f1=390;f2=410;wc1=2*f1/sf; wc2=2*f2/sf; %归一化频率f0=[0 wc1-0.05 wc1 wc2 wc2+0.05 1]; B=[0 0 1 1 0 0];%设置带通和带阻 weigh=[1 1 1 ];%设置带通和带阻权重 b=remez(50,f0,B,weigh);%传函分子 D=filter(b,1,N);y=sin(2*pi*t*100)+1.5*sin(2*pi*t*300)+A.*sin(2*pi*t*200)+D+N; a(1,:)=y;a(2,:)=y.*sin(y); x=a(1,:);y=a(2,:)-a(1,:);f=0:sf/nfft:sf/2-sf/nfft;w=boxcar(nfft);%加矩形窗 z=psd(y,nfft,sf,w,nfft/2); nn=1:nfft/2;plot(f(nn),abs(z(nn))); xlabel('频率(Hz)'); ylabel('幅值'); grid on;图4.1 功率谱估计结果图(2).信号频率的提取用离散傅立叶算法离散傅立叶算法程序 clear clcclose all hidden sf=1000;nfft=2048; t=0:1/1000:2.047;A=normrnd(0,1,1,2048); N=wgn(1,2048,1); f1=390;f2=410; wc1=2*f1/sf; wc2=2*f2/sf;050100150200250300350400450500200400600800频率(Hz)幅值%归一化频率f0=[0 wc1-0.05 wc1 wc2 wc2+0.05 1];B=[0 0 1 1 0 0];%设置带通和带阻weigh=[1 1 1 ];%设置带通和带阻权重b=remez(50,f0,B,weigh);%传函分子D=filter(b,1,N);y=sin(2*pi*t*100)+1.5*sin(2*pi*t*300)+A.*sin(2*pi*t*200)+D+N; t2=(0:nfft-1)/sf;f=(0:nfft-1)*sf/nfft;y1=abs(fft(y));f=f(1:nfft/2);y1=y1(1:nfft/2);plot(t,y);title('原始信号');axis([0 2.047 -6 8]);plot(f,y1);title('fft频率提取');axis([0 500 0 1000]);xlabel('f/Hz');grid on;原信号时间（t）图4.2 原始信号时域图图4.3 信号频谱(3)分别利用Winner 滤波和Kalman 滤波进行去噪；clear all close allM=100;%维纳滤波器阶数 sf=1000;nfft=2048; L=nfft;t=0:1/1000:2.047;A=normrnd(0,1,1,2048); N=wgn(1,2048,1); f1=390;f2=410; wc1=2*f1/sf; wc2=2*f2/sf; %归一化频率f0f0=[0 wc1-0.05 wc1 wc2 wc2+0.05 1]; B=[0 0 1 1 0 0];%设置带通和带阻 weigh=[1 1 1 ];%设置带通和带阻权重 b=remez(50,f0,B,weigh);%传函分子 D=filter(b,1,N);y=sin(2*pi*t*100)+1.5*sin(2*pi*t*300)+A.*sin(2*pi*t*200)+D+N; phixx=xcorr(y,y); for i=1:M for j=1:MRxx(i,j)=phixx(i-j+L); end ends=sin(2*pi*t*100)+1.5*sin(2*pi*t*300)+A.*sin(2*pi*t*200); phixs=xcorr(y,s); for i=1:Mrxs(i)=phixs(i+L); endh1=(inv(Rxx))*rxs';2004006008001000fft 频率提取f/Hz%获得理想FIR滤波器系数h1AA=sin(2*pi*t*100)+1.5*sin(2*pi*t*300)+A.*sin(2*pi*t*200); for i=1:Mh(i)=AA(i);end%绘图比较估计滤波器与实际滤波器figurek=1:M;plot(k,h(k),'r',k,h1(k),'b');title('Ideal h(n) & Calculated h(n)');legend('Ideal h(n)',' Calculated h(n)');xlabel('n');ylabel('h(n)');%比较理想输出与实际输出v=D+N;S=conv(h,v);SI(1)=S(1);LL1=sin(2*pi*t*100)+1.5*sin(2*pi*t*300)+A.*sin(2*pi*t*200); for i=2:LSI(i)=LL1(i);endfigurek=1:L;plot(k,s(k),'r',k,SI(k),'b');title('s(n) VS. SI(n)');legend('s(n)','SI(n)',0);xlabel('n');ylabel('Ideal Output');hold onSR=conv(h1,y);figurek=1:L;plot(k,s(k),'r',k,SR(k),'b');title('s(n)VS. SR(n)');legend('s(n)去噪前','SR(n)去噪后',0);xlabel('n');ylabel('Actual Output');图4.4 Winner 滤波去噪图Kalman 滤波程序 clear; clc;Fs=1000; nfft=2048;t1=0:1/Fs:2.047;A=normrnd(0,1,1,2048); N=wgn(1,2048,2); f1=390;f2=410; wc1=2*f1/Fs; wc2=2*f2/Fs; wc2=2*f2/sf; %归一化频率f0f0=[0 wc1-0.05 wc1 wc2 wc2+0.05 1]; B=[0 0 1 1 0 0];%设置带通和带阻 weigh=[1 1 1 ];%设置带通和带阻权重 b=remez(50,f0,B,weigh);%传函分子 D=filter(b,1,N);x=sin(2*pi*t1*100)+1.5*sin(2*pi*t1*300)+A.*sin(2*pi*t1*200)+D+N; x1=sin(2*pi*t1*100)+1.5*sin(2*pi*t1*300)+A.*sin(2*pi*t1*200); a1=-1.352;a2=1.338;a3=-0.662;a4=0.240;A=[-a1 -a2 -a3 -a4;1 0 0 0;0 1 0 0;0 0 1 0];%状态转移矩阵 H=[1 0 0 0];%观测矩阵Q=[1 0 0 0;0 0 0 0;0 0 0 0;0 0 0 0];%状态噪声方差 R=1;%观测噪声方差阵X(:,1)=[x(4);x(3);x(2);x(1)];p(:,:,1)=[10 0 0 0;0 1 0 0;0 0 1 0;0 0 0 1];%一步预测误差方针 %开始滤波 for k=2:nfftp1(:,:,k)=A*p(:,:,k-1)*A'+Q;%p1(:,:,k)即是一步预测误差的自相关矩阵，它是4*4的矩阵，取不同的k 值就构成了一个三维矩阵K(:,k)=p1(:,:,k)*H'/(H*p1(:,:,k)*H'+R); %K(:,:,k)是增益矩阵,对于固定的k 值它是4*1矩阵，取不同的k 值就是三维矩阵s(n)VS. SR(n)nA c t u a l O u t p u tX(:,k)=A*X(:,k-1)+K(:,k)*[x(k)-H*A*X(:,k-1)]; %X(:,k)是估计值,4*1矩阵p(:,:,k)=p1(:,:,k)-K(:,k)*H*p1(:,:,k);%p(:,:,k)是估计误差的自相关矩阵，4*4矩阵的三维矩阵end%结束一次滤波%绘图t=1:nfft;figure(2);plot(t,x1,'k-',t,x,'r-',t,X(1,:),'b-.');title('卡曼滤波去噪')legend('真实轨迹','观测样本','估计轨迹');grid on;卡曼滤波去噪n图5 Kalman滤波去噪图(4) 利用Wigner-Ville分布分析信号的时频特性MATLAB程序clear;clc;Fs=1000;nfft=2049;t1=0:1/Fs:2.048;A=normrnd(0,1,1,2049);N=wgn(1,2049,2);f1=390;f2=410;wc1=2*f1/Fs;wc2=2*f2/Fs;%归一化频率f0f0=[0 wc1-0.05 wc1 wc2 wc2+0.05 1];B=[0 0 1 1 0 0];%设置带通或带阻，1为带通，0为带阻weigh=[1 1 1 ];%设置通带和阻带的权重b=remez(50,f0,B,weigh);%传函分子D=filter(b,1,N);x=sin(2*pi*t1*100)+1.5*sin(2*pi*t1*300)+A.*sin(2*pi*t1*200)+D+N; figure(8)tfrwv(x');xlabel('时间t');ylabel('频率f');0.50.450.40.350.30.250.20.150.10.05图6 幅频特性图。
现代信号处理
4.信号的函数表达式为：()()()()sin(2100) 1.5sin(2300)sin(2200)x t t t A t t dn t n t πππ=++++，其中，()A t 为一随时间变化的随机过程，()dn t 为经过390—410Hz 带通滤波器后的高斯白噪声，()n t 为高斯白噪声，采样频率为1kHz ，采样时间为2.048s 。
（1）利用现代信号处理的知识进行信号谱估计； （2）利用现代信号处理知识进行信号的频率提取； （3）分别利用Winner 滤波和Kalman 滤波进行去噪； （4）利用Wigner-Ville 分布分析信号的时频特性。
（1）：利用现代信号处理的知识进行信号谱估计：经典谱估计中两种主要的方法为直接法和间接法，其中间接法则先根据N 个样本数据()x n 的样本自相关函数µ()()()1*01,01N x n Rk x n k x n k M N-==+=⋅⋅⋅∑，，，（4.1）其中1M N ≤<，且µ()µ()*x x R k R k -=。
计算样本自相关函数的Fourier 变换，得到功率谱()µ()Mjk x x k MP Rk e ωω-=-=∑（4.2）周期图方法估计的功率谱为有偏估计，可通过加窗来减少其偏差。
定义为 ()()()2101N jn x n P x n c n e NWωω--==∑ （4.3）式中()()122112N n W c n C d NNππωωπ--===∑⎰（4.4）式中，()C ω是窗函数()c n 的Fourier 变换。
功率谱估计程序为： clear clcclose all hidden sf=1000;nfft=2048; t=0:1/1000:2.047; A=normrnd(0,1,1,2048); N=wgn(1,2048,1); f1=390;f2=410; wc1=2*f1/sf; wc2=2*f2/sf; %归一化频率f0=[0 wc1-0.05 wc1 wc2 wc2+0.05 1]; B=[0 0 1 1 0 0];%设置带通和带阻 weigh=[1 1 1 ];%设置带通和带阻权重 b=remez(50,f0,B,weigh);%传函分子 D=filter(b,1,N);y=sin(2*pi*t*100)+1.5*sin(2*pi*t*300)+A.*sin(2*pi*t*200)+D+N; a(1,:)=y;a(2,:)=y.*sin(y); x=a(1,:); y=a(2,:)-a(1,:);f=0:sf/nfft:sf/2-sf/nfft; w=boxcar(nfft);%加矩形窗 z=psd(y,nfft,sf,w,nfft/2); nn=1:nfft/2;plot(f(nn),abs(z(nn))); xlabel('频率(Hz)'); ylabel('幅值'); grid on;图4.1 功率谱估计结果图(2).信号频率的提取用离散傅立叶算法离散傅立叶算法程序 clear clcclose all hidden sf=1000;nfft=2048; t=0:1/1000:2.047;050100150200250300350400450500200400600800频率(Hz)幅值A=normrnd(0,1,1,2048);N=wgn(1,2048,1);f1=390;f2=410;wc1=2*f1/sf;wc2=2*f2/sf;%归一化频率f0=[0 wc1-0.05 wc1 wc2 wc2+0.05 1];B=[0 0 1 1 0 0];%设置带通和带阻weigh=[1 1 1 ];%设置带通和带阻权重b=remez(50,f0,B,weigh);%传函分子D=filter(b,1,N);y=sin(2*pi*t*100)+1.5*sin(2*pi*t*300)+A.*sin(2*pi*t*200)+D+N; t2=(0:nfft-1)/sf;f=(0:nfft-1)*sf/nfft;y1=abs(fft(y));f=f(1:nfft/2);y1=y1(1:nfft/2);plot(t,y);title('原始信号');axis([0 2.047 -6 8]);plot(f,y1);title('fft频率提取');axis([0 500 0 1000]); xlabel('f/Hz'); grid on;图4.2 原始信号时域图图4.3 信号频谱(3)分别利用Winner 滤波和Kalman 滤波进行去噪；clear all close allM=100;%维纳滤波器阶数0.20.40.60.81 1.2 1.41.61.82原信号时间（t ）0501001502002503003504004505002004006008001000fft 频率提取f/Hzsf=1000;nfft=2048;L=nfft;t=0:1/1000:2.047;A=normrnd(0,1,1,2048);N=wgn(1,2048,1);f1=390;f2=410;wc1=2*f1/sf;wc2=2*f2/sf;%归一化频率f0f0=[0 wc1-0.05 wc1 wc2 wc2+0.05 1];B=[0 0 1 1 0 0];%设置带通和带阻weigh=[1 1 1 ];%设置带通和带阻权重b=remez(50,f0,B,weigh);%传函分子D=filter(b,1,N);y=sin(2*pi*t*100)+1.5*sin(2*pi*t*300)+A.*sin(2*pi*t*200)+D+N; phixx=xcorr(y,y);for i=1:Mfor j=1:MRxx(i,j)=phixx(i-j+L);endends=sin(2*pi*t*100)+1.5*sin(2*pi*t*300)+A.*sin(2*pi*t*200);phixs=xcorr(y,s);for i=1:Mrxs(i)=phixs(i+L);endh1=(inv(Rxx))*rxs';%获得理想FIR滤波器系数h1AA=sin(2*pi*t*100)+1.5*sin(2*pi*t*300)+A.*sin(2*pi*t*200); for i=1:Mh(i)=AA(i);end%绘图比较估计滤波器与实际滤波器figurek=1:M;plot(k,h(k),'r',k,h1(k),'b');title('Ideal h(n) & Calculated h(n)');legend('Ideal h(n)',' Calculated h(n)');xlabel('n');ylabel('h(n)');%比较理想输出与实际输出v=D+N;S=conv(h,v);SI(1)=S(1);LL1=sin(2*pi*t*100)+1.5*sin(2*pi*t*300)+A.*sin(2*pi*t*200);for i=2:LSI(i)=LL1(i);endfigurek=1:L;plot(k,s(k),'r',k,SI(k),'b');title('s(n) VS. SI(n)');legend('s(n)','SI(n)',0);xlabel('n');ylabel('Ideal Output'); hold onSR=conv(h1,y);figurek=1:L;plot(k,s(k),'r',k,SR(k),'b');title('s(n)VS. SR(n)');legend('s(n)去噪前','SR(n)去噪后',0); xlabel('n');ylabel('Actual Output');图4.4 Winner 滤波去噪图Kalman 滤波程序 clear; clc; Fs=1000; nfft=2048; t1=0:1/Fs:2.047; A=normrnd(0,1,1,2048); N=wgn(1,2048,2); f1=390;f2=410; wc1=2*f1/Fs; wc2=2*f2/Fs; wc2=2*f2/sf; %归一化频率f0f0=[0 wc1-0.05 wc1 wc2 wc2+0.05 1]; B=[0 0 1 1 0 0];%设置带通和带阻 weigh=[1 1 1 ];%设置带通和带阻权重s(n)VS. SR(n)nA c t u a l O u t p u tb=remez(50,f0,B,weigh);%传函分子D=filter(b,1,N);x=sin(2*pi*t1*100)+1.5*sin(2*pi*t1*300)+A.*sin(2*pi*t1*200)+D+N;x1=sin(2*pi*t1*100)+1.5*sin(2*pi*t1*300)+A.*sin(2*pi*t1*200);a1=-1.352;a2=1.338;a3=-0.662;a4=0.240;A=[-a1 -a2 -a3 -a4;1 0 0 0;0 1 0 0;0 0 1 0];%状态转移矩阵H=[1 0 0 0];%观测矩阵Q=[1 0 0 0;0 0 0 0;0 0 0 0;0 0 0 0];%状态噪声方差R=1;%观测噪声方差阵X(:,1)=[x(4);x(3);x(2);x(1)];p(:,:,1)=[10 0 0 0;0 1 0 0;0 0 1 0;0 0 0 1];%一步预测误差方针%开始滤波for k=2:nfftp1(:,:,k)=A*p(:,:,k-1)*A'+Q;%p1(:,:,k)即是一步预测误差的自相关矩阵，它是4*4的矩阵，取不同的k值就构成了一个三维矩阵K(:,k)=p1(:,:,k)*H'/(H*p1(:,:,k)*H'+R); %K(:,:,k)是增益矩阵,对于固定的k 值它是4*1矩阵，取不同的k值就是三维矩阵X(:,k)=A*X(:,k-1)+K(:,k)*[x(k)-H*A*X(:,k-1)]; %X(:,k)是估计值,4*1矩阵p(:,:,k)=p1(:,:,k)-K(:,k)*H*p1(:,:,k);%p(:,:,k)是估计误差的自相关矩阵，4*4矩阵的三维矩阵end%结束一次滤波%绘图t=1:nfft;figure(2);plot(t,x1,'k-',t,x,'r-',t,X(1,:),'b-.');title('卡曼滤波去噪')legend('真实轨迹','观测样本','估计轨迹');grid on;卡曼滤波去噪n图5 Kalman滤波去噪图(4) 利用Wigner-Ville分布分析信号的时频特性MATLAB程序clear;clc;Fs=1000;nfft=2049;t1=0:1/Fs:2.048;A=normrnd(0,1,1,2049);N=wgn(1,2049,2);f1=390;f2=410;wc1=2*f1/Fs;wc2=2*f2/Fs;%归一化频率f0f0=[0 wc1-0.05 wc1 wc2 wc2+0.05 1];B=[0 0 1 1 0 0];%设置带通或带阻，1为带通，0为带阻weigh=[1 1 1 ];%设置通带和阻带的权重b=remez(50,f0,B,weigh);%传函分子D=filter(b,1,N);x=sin(2*pi*t1*100)+1.5*sin(2*pi*t1*300)+A.*sin(2*pi*t1*200)+D+N; figure(8)tfrwv(x');xlabel('时间t');ylabel('频率f');0.50.450.40.350.30.250.20.150.10.05图6 幅频特性图。
机械故障诊断中的现代信号处理方法
机械故障诊断中的现代信号处理方法
现代信号处理方法在机械故障诊断中有着广泛的应用。
以下是几种常见的现代信号处理方法：
1. 傅里叶变换（Fourier Transform）: 傅里叶变换将时域信号转换为频域信号，可以分析信号的频率成分和能量分布。
在机械故障诊断中，傅里叶变换可以用来检测故障产生的谐波或频率成分的变化。
2. 小波变换（Wavelet Transform）: 小波变换可以在时间和频率上同时进行分析，可以更好地捕捉瞬态故障或频率变化的特征。
小波变换在机械故障诊断中常用于检测冲击、噪声和频率模态等问题。
3. 自适应滤波（Adaptive Filtering）: 自适应滤波是一种可以自动调整滤波器参数的方法，可以根据信号的特点动态调整滤波器的频率响应。
自适应滤波在机械故障诊断中可以用于降噪和提取故障特征。
4. 统计特征提取（Statistical Feature Extraction）: 统计特征提取是通过对信号进行统计分析来提取信号特征的方法。
常见的统计特征包括均值、方差、峰值、峭度等。
统计特征提取可以用来检测信号的变化和异常。
5. 机器学习（Machine Learning）: 机器学习是一种可以让计算机自动学习和适应数据模式的方法。
在机械故障诊断中，机器学习可以用来训练模型，识别和分类不同的故障模式。
常见的
机器学习算法包括支持向量机（SVM）、随机森林（Random Forest）和深度学习（Deep Learning）等。
这些现代信号处理方法可以结合使用，以提取和分析机械故障信号中的相关特征，提高故障诊断的准确性和效率。
现代信号处理技术
与功率谱分析比较，运用基于高阶累计量的谱估 计算法估计信号，消除了高斯噪声的影响，使估计结 果更准确，并且保留了信号的相位特性，提供更多的 内在信息。
四、Choi-Williams 分布(CWD)
WD分布来源于广义时频分布，定义为：
C W D (t, )
e (4 t 2 u )2x (u)x * (u)e j d u d
42
22
(11-3)
通常,在处理幅度和频率变化较大的信号时取较大的R(R>1) 值;反之,则取较小R(R≤1) 值。CWD满足多数所希望的时 频特性,其抑制交叉项的能力还取决于被分析信号的时频 结构。因此,实际应用中需要综合考虑。
五、Cone 核分布(CKD) 等
当核函数 (t,) 1e2
0 一步变成Cone核分布：
t 时，广义时频分布进
其它
C K D (t, ) 1 e 2x(u)x*(u)ej d u d
22
(11-4)
式中， t 。
CKD 具有较好的抑制横向交叉项的能力, 适合处理这样的 信号, 即在一个小的范围内频率分布是正值, 而在此之外频 率分布是负值, 参数R确定范围的大小。
得模糊。
三、Wigner-Ville 分布(WVD)
实际信号s(t) 的Wigner-Ville 分布定义为:
W V D (t,) x(t)x*(t)ej d
2 2
(11-2)
式中: x(t)为s(t)的解析信号。 在Wigner-Ville 分布中使用解析信号x(t)而不是 原实际信号s(t)的优点在于: 第一,解析信号的处 理中只采用频谱正半部分,因此不存在由正频率 项和负频率项产生的交叉项；第二,使用解析信 号不需要过采样,同时可避免不必要的畸变影响。
现代信号处理
2015年12月20日
机械工程学院机自所动态室
3
第七章 基于第二代小波变换的信号处理
7.1 第二代小波变换原理 7.2 预测器和更新器 7.3 第二代小波包分析 7.4 冗余第二代小波变换
2015年12月20日
机械工程学院机自所动态室
4
第七章 基于第二代小波变换的信号处理
7.1 第二代小波变换原理 7.2 预测器和更新器 7.3 第二代小波包分析 7.4 冗余第二代小波变换
(7.1.5)
右边界受影响的情况有 D 种，预测器统一表示为
P(se ) p1se (L' N 1) p2se (L' N 2) ... pN se (L' )
(7.1.6)
L'为偶样本序列 se 的长度。
2015年12月20日
机械工程学院机自所动态室
8
7.1 第二代小波变换原理
(7.2.3)
2015年12月20日
机械工程学院机自所动态室
16
7.2 预测器和更新器
7.2.2 更新器系数计算方法
设在更新阶段，更新器U 的个数为 N~（N~ 2D~ ，D~ 为正整 数），预测器 P 的个数为N（N 2D ，D 为正整数）。
将 P 和 U 代入第二代小波重构等效高通滤波器表达式， 则得到重构等效高通滤波器 g表达式如下
2015年12月20日
机械工程学院机自所动态室
10
7.1 第二代小波变换原理
第二代小波变换的重构过程由三部分组成：恢复更新、
恢复预测和合并。其过程实现如图7.1.3所示。
s
se
重构
-U
P
(merge)
S
d
现代信号处理
课程简介
现代信号处理是“信息与通信工程 信息与通信工程”一级学 科“通信与信息系统”和“信号与信息处理” 两个专业的学位课，“电子与通信工程 电子与通信工程”专 业 门重要的学位专业基 课 作为信息处 业一门重要的学位专业基础课，作为信息处 理与现代通信的基础，它对信息科学的发展 起着重要作用 先修课程：随机过程、最优化方法、数字信 先修课程 号处理
现代信号处理 8
lifei@
信号处理的典型应用
• • • • • • 1.语音处理 2.图像处理 3.通信 4 航空航天 4.航空航天 5.生物医学 ……
语音处理
• • • • • 最早采用DSP的领域之一 语音编码 语音合成 语音识别 语音增强
lifei@
lifei@
图像处理
• 数据压缩 • 图像复原 • 清晰化与增强
信号处理方法（小结）
• 方法分类
– – – – 基于变换的方法（Fourier 变换） 基于统计的方法 (Bayes准则) 基于模型的方法 (信号模型AR, AR MA MA, ARMA) 基于智能/机器学习的方法 (盲方法，对信号所知甚少)
现代信号处理 - 研究内容
DSP DSP: 两大支柱,表层信息 –快速变换 –数字滤波 MSP MSP: 四大处理, 深层信息 –自适应信号处理(盲,半盲） –现代谱估计（如HOS） –非平稳信号处理(Wavelets) –非线性信号处理(如NNSP)
现代信号处理 20
lifei@
现代信号处理 - 处理方法
• 取决于信号本身（关于对信号本身的知识） • 取决于具体应用
信号处理框图
D S P
现代信号处理 21
ห้องสมุดไป่ตู้
现代信号处理方法自适应信号处理方法
yj XT jWWTXj
式中
(2.1.3)
W [ w 1 ,w 2 , ,w N ] T ,X j [ x 1 j,x 2 j, ,x N ] T j
误差信号表示为 e j d j y j d j X T jW d j W T X j
(2.1.4)
现代信号处理方法自适应信号处理 方法
自适应信号处理
现代信号处理方法自适应信号处理 方法
自适应信号处理
二、LMS 自适应横向滤波器 LMS自适应滤波器是以均方误差最小作为最佳滤波
准则的，原理框图如图2.1所示，图中x(n)称为输入信 号，y(n)是输出信号，d(n)称为期望信号，或者称为参
考信号、训练信号，e(n)是误差信号。
e(n)=d(n)-y(n)
x(n) H(z)
e(n)
y(n)
－
＋
d(n)
现代信号处图理方2法.1自适应自信号适处应理 滤波器原理图
方法
自适应信号处理
二、LMS 自适应横向滤波器 其中自适应滤波器H(z)的系数根据误差信号，通过
一定的自适应算法，不断地进行改变,使输出y(n)最接近 期望信号d(n),这里暂时假定d(n)是可以利用的，实际中， d(n)要根据具体情况进行选取,能够选到一个合适的信 号作为期望信号，是设计自适应滤波器的一项重要的 工作。如果真正的d(n)可以获得， 我们将不需要做任 何自适应滤波器。
…
x1j
w1
x2j
w2
xNj wN
yj
－
＋
ej
dj
图 2 现自代适信号应处理线方性法自组适应合信号器处理
方法
自适应信号处理
二、LMS 自适应横向滤波器 2.1. 自适应滤波器的矩阵表示式
专业学位硕士研究生“现代信号处理”课程教学改革探讨
专业学位硕士研究生“现代信号处理”课程教学改革探讨一、课程教学内容和目标“现代信号处理”是一门与电子信息工程相关的专业课程，旨在使学生掌握信号处理的基础知识和技术方法，能够熟练掌握信号处理的基本理论和方法，具备基本的信号处理系统设计和实现能力。
而在教学内容方面，应该包括信号与系统的基本概念、离散信号的表示与处理、傅里叶变换与频谱分析、滤波器设计与实现等内容。
目前，传统的信号处理课程教学过于理论化，教学内容较为单一，无法满足学生的实际需求。
需要对“现代信号处理”课程的教学内容和目标进行重新规划和调整，以适应时代的发展和改变。
在今后的课程教学改革中，应该重点加强对信号处理的实际应用和工程技术方面的教学，使学生能够更好地掌握信号处理的基础理论和应用技能，提高学生的实践能力和创新思维。
二、教学方法和手段教学方法和手段是课程教学改革中需要着重考虑和调整的方面。
传统的课堂教学方式已经无法满足学生的学习需求，我们需要引入更多的实践教学和案例教学，通过项目实践、实验教学等方式，来培养学生的创新能力和实际操作能力。
在进行课程教学改革时，可以适当增加一些与实际工程应用相关的案例分析和真实工程项目，以提高学生的学习兴趣，并激发学生的学习激情。
在教学手段方面，可以适当引入多媒体教学、网络教学等现代化教学手段，以提高教学效果和学生的学习兴趣。
通过使用多媒体教学和网络教学工具，可以使学生更直观地了解现代信号处理技术的应用，增强学生的实践能力和创新意识。
三、课程教学评估课程教学的评估是课程教学改革工作的关键环节，对教学改革的成效进行及时和科学的评估，对今后的教学改革工作起着十分重要的作用。
在进行课程教学评估时，应该从多个方面进行评价，包括学生的学习效果、教学方法的有效性、教学手段的使用效果等。
在对学生的学习效果进行评估时，可以采取考试成绩、课程作业、实习报告等多种评价手段，全面地评价学生在“现代信号处理”课程学习过程中的表现。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
一、有两个ARMA 过程，其中信号1是宽带信号，信号2是窄带信号，分别用AR 谱估计算法、ARMA 谱估计算法和周期图算法估计其功率谱。
产生信号1的系统函数为1234123410.35440.35080.17360.2401()1 1.3817 1.56320.88430.4906z z z z H z z z z z --------++++=-+-+激励白噪声的方差为1。
产生信号2的系统函数为1212341 1.58570.9604()1 1.6408 2.2044 1.48080.8145z z H z z z z z ------++=-+-+激励白噪声的方差为1。
每次实验使用的数据长度为256。
(1) 对信号1，分别使用AR(4)，AR(8)，ARMA(4,4)和ARMA(8,8)模型进行谱估计，对AR 方法采用自协方差算法，对ARMA 算法采用改进的Yule-Walker 方程算法，也用周期图法作谱估计。
做20次独立实验，将20次实验结果画在一张图上，观察谱估计的随机分布性质，另将20次的平均值和真实谱画在一张图上进行比较。
(2) 对信号2，分别使用AR(4)，AR(8)，AR(12)，AR(16)，ARMA(4,2)，ARMA(8,4)和ARMA(12,6)模型进行谱估计，对AR 方法采用自协方差算法，对ARMA 算法采用改进的Yule-Walker 方程算法，也用周期图法作谱估计。
做20次独立实验，将20次实验结果画在一张图上，观察谱估计的随机分布性质，另将20次的平均值和真实谱画在一张图上进行比较。
(3) 对各种算法的性能进行比较分析。
解:(1 )很多随机过程可以由或近似由均值为零、方差为 的白噪声序列u(n)经过具有有理传输函数H(z)的ARMA 线性系统来得到。
称该随机过程为ARMA 过程。
任意平稳ARMA 过程，其功率谱密度有如下形式：222|)(||)(|)(σωωω⋅=A B P x (1)则x(n)可用ARMA(p,q)模型描述，即22|)(|)(σωω⋅=H P x(2)则可以根据给出的信号1的系统函数来进行计算。
2σ仿真图形结果下所示：信号1的20次谱估计总图如图一所示：信号1的20次AR(4)谱估计总图A B信号1的20次的ARMA(4,4)谱估计总图信号1的20次的ARMA(8,8)谱估计总图C D 图一 信号1的20次谱估计信号1的20次谱估计平均如下图二所示：信号1的20次AR4谱估计平均图信号1的20次AR8谱估计平均a b信号1的20次ARMA(4,4)谱估计平均图信号1的20次ARMA(8,8)谱估计平均C d图二信号1的20次谱估计平均上图中蓝色曲线表示理论值，红色曲线表示估计值。
信号1的20次周期图法图谱估计如图三所示：信号1的20次周期图法谱估计平均图a b图三信号1的20次周期图法图(2)仿真结果如下所示：信号2的20次谱估计总图如图四所示信号2的20次AR(4)谱估计总图信号2的20次AR(8)谱估计总图a bC de fG图四信号2的20次谱估计信号2的20次谱估计平均如图五所示：信号2的20次AR4谱估计平均信号2的20次AR8谱估计平均信号2的20次AR12谱估计平均信号2的20次ARMA(4,2)谱估计平均信号2的20次ARMA(8,4)谱估计平均信号2的20次ARMA(12,6)谱估计平均图六 信号2的20次谱估计平均信号2的周期法谱估计如图七所示：信号2的20次的周期图法谱估计总图信号2的20次周期图法谱估计平均3 各个算法分析，对于信号1，有图形可以看出，AR模型，估计所得值随着阶数的增加，更加靠近真实值。
ARMA算法，阶数的增加，估计曲线与理论曲线更加吻合。
相同条件下ARMA比AR所得的效果更好。
从20次平均结果来看，周期图法最接近理论值。
对于信号2得到的结果和信号1类似。
MATLAB程序代码信号1clear all;close all;clc;N = 256;omiga = pi/N:pi/N:pi;A = [1 -1.3817 1.5632 -0.8843 .4906];B = [1 .3544 .3508 .1736 .2401];Px_peri = zeros(20,N);Px_peritotal = zeros(1,N);Px_ar4 = zeros(20,N);Px_ar4total = zeros(1,N);Px_arma44 = zeros(20,N);Px_arma44total = zeros(1,N);Px_ar8 = zeros(20,N);Px_ar8total = zeros(1,N);Px_arma88 = zeros(20,N);Px_arma88total = zeros(1,N);Px_arth = get_pw(A(2:5),0,4,1,N);Px_armath = get_pw(A(2:5),B(2:5),4,4,N);for l = 1:20u = normrnd(0,1,1,N);x = filter(B,A,u);p = 4;A_mation4 = get_ar_coe(p,x);Px_ar4(l,:) = get_pw(A_mation4,0,p,1,N);q = 4;y = filter([1 A_mation4],1,x);p_temp = 100;a_temp = get_ar_coe(p_temp,y);B_mation4 =get_ar_coe(q,a_temp);Px_arma44(l,:) = get_pw(A_mation4,B_mation4,p,q,N);p = 8;A_mation8 = get_ar_coe(p,x);Px_ar8(l,:) = get_pw(A_mation8,0,p,1,N);q = 8;y = filter([1 A_mation8],1,x);p_temp = 100;a_temp = get_ar_coe(p_temp,y);B_mation8 =get_ar_coe(q,a_temp);Px_arma88(l,:) = get_pw(A_mation8,B_mation8,p,q,N);for j = 1:NPx_peri(l,j) = sum(x.*exp(-1i*omiga(j)*(1:256)))*conj(sum(x.*exp(-1i*omiga(j)*(1:256))));endendPx_peri = Px_peri/N;figure;hold on;for l = 1:20plot(omiga,log(Px_ar4(l,:)));Px_ar4total = Px_ar4total + Px_ar4(l,:);endPx_ar4total = Px_ar4total/20;title('信号1的20次AR(4)谱估计总图');axis([0 pi -5 5]);figure;plot(omiga,log(Px_armath));title('信号1的ARMA模型理论谱');axis([0 pi -5 5]);hold on;plot(omiga,log(Px_ar4total),'-r');title('信号1的20次AR4谱估计平均图');axis([0 pi -5 5]);figure;hold on;for l = 1:20plot(omiga,log(Px_ar8(l,:)));Px_ar8total = Px_ar8total + Px_ar8(l,:);endPx_ar8total = Px_ar8total/20;title('20次信号1的AR(8)谱估计总图');axis([0 pi -5 5]);figure;plot(omiga,log(Px_armath));title('信号1的ARMA模型理论谱');axis([0 pi -5 5]);hold on;plot(omiga,log(Px_ar8total),'-r');title('信号1的20次AR8谱估计平均');axis([0 pi -5 5]);figure;hold on;for l = 1:20plot(omiga,log(Px_arma44(l,:)));Px_arma44total = Px_arma44total + Px_arma44(l,:); endPx_arma44total = Px_arma44total/20;title('信号1的20次的ARMA(4,4)谱估计总图');axis([0 pi -5 5]);figure;plot(omiga,log(Px_armath));title('信号1的ARMA模型理论谱图');axis([0 pi -5 5]);hold on;plot(omiga,log(Px_arma44total),'-r');title('信号1的20次ARMA(4,4)谱估计平均图');axis([0 pi -5 5]);figure;hold on;for l = 1:20plot(omiga,log(Px_arma88(l,:)));Px_arma88total = Px_arma88total + Px_arma88(l,:); endPx_arma88total = Px_arma88total/20;title('信号1的20次的ARMA(8,8)谱估计总图');axis([0 pi -5 5]);figure;plot(omiga,log(Px_armath));title('信号1的ARMA模型理论谱');axis([0 pi -5 5]);hold on;plot(omiga,log(Px_arma88total),'-r');title('信号1的20次ARMA(8,8)谱估计平均');axis([0 pi -5 5]);figure;hold on;for l = 1:20plot(omiga,log(Px_peri(l,:)));Px_peritotal = Px_peritotal + Px_peri(l,:);endPx_peritotal = Px_peritotal/20;title('信号1的20次的周期图法谱估计总图');axis([0 pi -5 5]);figure;plot(omiga,log(Px_armath));title('信号1的ARMA模型理论谱图');axis([0 pi -5 5]);hold on;plot(omiga,log(Px_peritotal),'-r');title('信号1的20次周期图法谱估计平均图');axis([0 pi -5 5]);信号2clear all;close all;clc;N = 256; % 信号点数omiga = pi/N:pi/N:pi; % 谱估计x轴坐标B = [1 1.5857 .9604];A = [1 -1.6408 2.2044 -1.4808 .8145];Px_peri = zeros(20,N);Px_peritotal = zeros(1,N);Px_ar4 = zeros(20,N);Px_ar4total = zeros(1,N);Px_arma42 = zeros(20,N);Px_arma42total = zeros(1,N); Px_ar8 = zeros(20,N);Px_ar8total = zeros(1,N);Px_arma84 = zeros(20,N);Px_arma84total = zeros(1,N); Px_ar12 = zeros(20,N);Px_ar12total = zeros(1,N);Px_arma126 = zeros(20,N);Px_arma126total = zeros(1,N); Px_ar16 = zeros(20,N);Px_ar16total = zeros(1,N);Px_arth = get_pw(A(2:5),0,4,1,N);Px_armath = get_pw(A(2:5),B(2:3),4,2,N);for l = 1:20u = normrnd(0,1,1,N);x1 = filter(B,A,u);p = 4;A_mation4 = get_ar_coe(p,x1);Px_ar4(l,:) = get_pw(A_mation4,0,p,1,N);q = 2;y1 = filter([1 A_mation4],1,x1);p_temp = 100;a_temp = get_ar_coe(p_temp,y1);B_mation2 =get_ar_coe(q,a_temp);Px_arma42(l,:) = get_pw(A_mation4,B_mation2,p,q,N);p = 8;A_mation8 = get_ar_coe(p,x1);Px_ar8(l,:) = get_pw(A_mation8,0,p,1,N);q = 4;y1 = filter([1 A_mation8],1,x1);p_temp = 100;a_temp = get_ar_coe(p_temp,y1);B_mation4 =get_ar_coe(q,a_temp);Px_arma84(l,:) = get_pw(A_mation8,B_mation4,p,q,N);p = 12;A_mation12 = get_ar_coe(p,x1);Px_ar12(l,:) = get_pw(A_mation12,0,p,1,N);q = 6;y1 = filter([1 A_mation12],1,x1);p_temp = 100; % 取阶数为100a_temp = get_ar_coe(p_temp,y1); %求Y1(z)等价100阶AR模型的解B_mation6 =get_ar_coe(q,a_temp); % 由a_temp建立q=6阶AR模型Px_arma126(l,:) = get_pw(A_mation12,B_mation6,p,q,N);% 得到ARMA(12,6)的谱估计p = 16;A_mation16 = get_ar_coe(p,x1); % 得到AR(16)的估计参数Px_ar16(l,:) = get_pw(A_mation16,0,p,1,N); % 得到AR(16)的谱估计for j = 1:NPx_peri(l,j) = sum(x1.*exp(-1i*omiga(j)*(1:256)))*conj(sum(x1.*exp(-1i*omiga(j)*(1:256))));endendPx_peri = Px_peri/N;figure;hold on;for l = 1:20plot(omiga,log(Px_ar4(l,:)));Px_ar4total = Px_ar4total + Px_ar4(l,:); endPx_ar4total = Px_ar4total/20;title('信号2的20次AR(4)谱估计总图');axis([0 pi -12 10]);figure;plot(omiga,log(Px_armath));title('信号2的ARMA模型理论谱');axis([0 pi -12 10]);hold on;plot(omiga,log(Px_ar4total),'-r');title('信号2的20次AR4谱估计平均');axis([0 pi -12 10]);figure;hold on;for l = 1:20plot(omiga,log(Px_ar8(l,:)));Px_ar8total = Px_ar8total + Px_ar4(l,:); endPx_ar8total = Px_ar8total/20;title('信号2的20次AR(8)谱估计总图');axis([0 pi -12 10]);figure;plot(omiga,log(Px_armath));title('信号2的ARMA模型理论谱');axis([0 pi -12 10]);hold on;plot(omiga,log(Px_ar8total),'-r');title('信号2的20次AR8谱估计平均');axis([0 pi -12 10]);figure;hold on;for l = 1:20plot(omiga,log(Px_ar12(l,:)));Px_ar12total = Px_ar12total + Px_ar4(l,:); endPx_ar12total = Px_ar12total/20;title('信号2的20次AR(12)谱估计总图');axis([0 pi -12 10]);figure;plot(omiga,log(Px_armath));title('信号2的ARMA模型理论谱');axis([0 pi -12 10]);hold on;plot(omiga,log(Px_ar12total),'-r');title('信号2的20次AR12谱估计平均');axis([0 pi -12 10]);figure;hold on;for l = 1:20plot(omiga,log(Px_ar16(l,:)));Px_ar16total = Px_ar16total + Px_ar4(l,:);endPx_ar16total = Px_ar16total/20;title('信号2的20次AR(16)谱估计总图');axis([0 pi -12 10]);figure;plot(omiga,log(Px_armath));title('信号2的ARMA模型理论谱');axis([0 pi -12 10]);hold on;plot(omiga,log(Px_ar16total),'-r');title('信号2的20次AR16谱估计平均');axis([0 pi -12 10]);figure;hold on;for l = 1:20plot(omiga,log(Px_arma42(l,:)));Px_arma42total = Px_arma42total + Px_arma42(l,:); endPx_arma42total = Px_arma42total/20;title('信号2的20次的ARMA(4,2)谱估计总图');axis([0 pi -12 10]);figure;plot(omiga,log(Px_armath));title('信号2的ARMA模型理论谱');axis([0 pi -12 10]);hold on;plot(omiga,log(Px_arma42total),'-r');title('信号2的20次ARMA(4,2)谱估计平均');axis([0 pi -12 10]);figure;hold on;for l = 1:20plot(omiga,log(Px_arma84(l,:)));Px_arma84total = Px_arma84total + Px_arma84(l,:); endPx_arma84total = Px_arma84total/20;title('信号2的20次的ARMA(8,4)谱估计总图');axis([0 pi -12 10]);figure;plot(omiga,log(Px_armath));title('信号2的ARMA模型理论谱');axis([0 pi -12 10]);hold on;plot(omiga,log(Px_arma84total),'-r');title('信号2的20次ARMA(8,4)谱估计平均');axis([0 pi -12 10]);figure;hold on;for l = 1:20plot(omiga,log(Px_arma126(l,:)));Px_arma126total = Px_arma126total + Px_arma126(l,:); endPx_arma126total = Px_arma126total/20;title('信号2的20次的ARMA(12,6)谱估计总图');axis([0 pi -12 10]);figure;plot(omiga,log(Px_armath));title('信号2的ARMA模型理论谱');axis([0 pi -12 10]);hold on;plot(omiga,log(Px_arma126total),'-r');title('信号2的20次ARMA(12,6)谱估计平均');axis([0 pi -12 10]);figure;hold on;for l = 1:20plot(omiga,log(Px_peri(l,:)));Px_peritotal = Px_peritotal + Px_peri(l,:); endPx_peritotal = Px_peritotal/20;title('信号2的20次的周期图法谱估计总图');axis([0 pi -12 10]);figure;plot(omiga,log(Px_armath));title('信号2的ARMA模型理论谱');axis([0 pi -12 10]);hold on;plot(omiga,log(Px_peritotal),'-r');title('信号2的20次周期图法谱估计平均');axis([0 pi -12 10]);。